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Recent  techniques  in  nonlinear  vibro-acoustics  have  demonstrated  improved  sens¬ 
ing  capabilities  for  landmine  detection.  These  methods  however,  place  the  transmit  and/or 
receive  devices  extremely  close  to  a  potentially  dangerous  target.  This  paper  discusses  a 
novel  approach  where  ultrasonic  parametric  arrays  are  used  to  achieve  excitation  at  standoff 
ranges  in  air.  When  two  frequencies,  f\  and  f-2  are  directed  to  excite  a  target,  the  nonlinear 
response  consists  of  sum  and  difference  frequencies. 

The  difference  frequency  may  be  carefully  swept  to  produce  an  acoustic  signature 
of  the  target,  reflecting  its  size  and  density  information.  As  part  of  this  research,  a  more 
accurate  third  order  nonlinear  ultrasonic  propagation  model  is  developed  to  analyze  signal 
strength  and  frequency  at  the  target.  Due  to  the  inefficient  mixing  of  the  ultrasonic  tones, 
reflected  signals  have  very  small  amplitude.  This  thesis  develops  high-resolution  spectral 
analysis  techniques  (e.g.  multiple  signal  classification  (MUSIC)  algorithm)  to  extract  par¬ 
ticularly  weak  signals  (in  low  signal  to  noise  ratio  scenario)  and  thus  substantially  improve 
target  characteristics  estimation  performance  and  provides  a  viable  and  practical  approach 
to  perform  acoustic  imaging.  Experimental  results  demonstrate  for  the  first  time,  a  capacity 
to  remotely  classify  a  hollow  target  from  a  solid  one,  with  resonance  patterns  predicting  the 
approximate  size  of  the  target. 
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Chapter  1 

Introduction 

1.1  Motivation 

The  use  of  sound,  especially  at  ultrasonic  frequencies,  has  been  used  extensively  for 
the  purposes  of  biological  imaging  and  non-destructive  testing.  These  methods  are  limited 
in  their  imaging  abilities  for  two  key  reasons:  in  non-destructive  testing,  the  transducer  must 
be  coupled  to  the  test  material  by  a  liquid  medium  for  impedance  matching  so  as  to  transmit 
as  much  energy  as  possible  into  the  test  object.  In  addition  to  this  limitation,  ultrasound 
is  only  able  to  measure  range  by  calculating  the  time  delay  of  a  reflected  signal.  Such 
techniques  rely  on  orthogonal  reflections  from  acoustically  reflective  surfaces.  Furthermore, 
there  is  inherent  noise  in  linear  acoustic  systems  because  the  reflected,  measured  acoustic 
energy  is  often  masked  by  the  transmitted  energy,  which  is  always  at  the  same  frequency. 

Recently,  techniques  that  use  two  frequencies,  and  the  nonlinear  effect  in  air  or 
other  media  have  shown  promise  in  gathering  even  more  information  about  our  environment. 
This  information  may  include  not  only  range  and  shape  data,  but  also  resonance  and  density 
measurements.  This  can  lead  to  better  examination  techniques  for  buried  landmines  and 
archeological  artifacts,  which  form  the  focus  of  this  thesis.  Furthermore,  research  in  this 
field  can  improve  techniques  for  air  coupled  ultrasonic  inspection  where  a  large  impedance 
mismatch  limits  the  ability  to  transmit  high-energy  signals  from  one  media  to  another  [5]. 

Current  technologies  that  are  used  for  the  detection  of  buried  objects  rely  on  the 
transmission  and  measurement  of  directed  energy.  These  methods  include  ground  penetrat¬ 
ing  radar,  infrared,  neutron  activation  analysis,  and  acoustics[6].  Most  of  these  technologies, 
however,  are  limited  in  their  ability  to  differentiate  targets  from  other  debris  in  the  soil. 
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Furthermore,  they  must  heavily  rely  on  assumptions  about  the  target  or  surrounding  media 
such  as  density,  compliance,  and  moisture.  Donskoy[7]  proposed  a  nonlinear  acoustic  tech¬ 
nique  for  landmine  detection  in  which  two  tones  are  used  to  insonify  the  soil.  Here  insonify 
is  used  to  describe  the  process  of  exposing  the  target  to  acoustic  or  sonar  energy  [8].  The 
key  assumption  to  this  approach  is  that  the  top  plate  of  the  landmine  must  have  a  stiffness 
less  than  or  equal  to  that  of  the  surrounding  soil.  Because  of  the  compliant  nature  of  the 
landmine  top  plate,  a  nonlinear  interaction  is  created  between  it  and  the  soil  above.  It  is  this 
nonlinear  interaction,  shown  in  Figure  1.1,  which  creates  a  difference  frequency,  indicating 
the  presence  of  a  buried  compliant  object  or  landmine.  This  technique  is  superior  to  other 


3.5"  Plastic  Mine,  1 "  depth 


Figure  1.1:  Acoustic  spectrum  of  3.5”  plastic  landmine.  Excitation  signals  are  visible  at  f\ 
=  400  Hz  and  f-2  =  650  Hz.  The  nonlinearity  at  f-2  —  fi  =  250  Hz  indicates  the  presence  of 
a  landmine  [1]. 

current  forms  of  acoustic  detection  because  the  excitation  signals  are  not  needed  in  the  de¬ 
tection  stage.  In  linear  detection,  information  about  the  target  is  contained  in  the  reflection 
of  the  excitation  signal.  This  reflected  energy  is  usually  very  weak  and  easily  masked  by 
the  incident  energy  as  well  as  reflections  from  many  non-target  objects  such  as  the  surface 
of  the  soil,  rocks,  inhonrogeneities,  and  other  battlefield  debris.  Because  the  reflected  signal 
is  at  a  different  frequency  than  the  excitation  signals,  nonlinear  vibro-acoustic  techniques 
are  much  more  discriminating  than  other  techniques.  Although  Donskoy  and  Kornran[2] 
have  demonstrated  the  viability  of  this  technique,  both  rely  on  geophones  and  loudspeakers 
that  must  be  placed  dangerously  close  to  the  target.  Figure  1.2  shows  Kornran’s  test  setup 
where  the  loudspeakers  and  geophones  must  be  placed  directly  above  the  landmine. 
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Figure  1.2:  Test  setup  for  Korman’s  2002  nonlinear  two-tone  acoustic  test  for  VS  2.2  land¬ 
mine.  The  mine  is  buried  at  a  depth  of  3.6  cm  and  located  7.9  cm  above  the  concrete  base. 
The  loudspeakers,  microphone,  and  geophone  must  be  placed  directly  above  the  mine[2]. 

1.2  Contribution 

This  thesis  proposes  the  use  of  highly  directional  ultrasonic  parametric  arrays  to 
provide  nonlinear  excitation  of  targets  at  standoff  ranges.  Using  the  concept  of  “scattering  of 
sound  by  sound”,  we  can  use  two  pure  tones,  /i  and  /2,  to  generate  a  difference  frequency  in 
air.  This  concept  was  presented  by  Yoneyama[4]  and  is  discussed  in  greater  detail  in  Section 
1.3.4.  This  difference  frequency  may  be  swept  to  provide  a  frequency  response  or  acoustic 
signature  of  a  target.  The  acoustic  signature  indicates  at  which  frequencies  a  target  absorbs 
and  reflects  energy  and  possibly  the  existence  of  resonant  behaviors.  While  detecting  buried 
landmines  is  still  beyond  the  scope  of  this  research,  discriminating  between  solid  and  hollow 
targets,  and  observing  resonant  behaviors  can  still  provide  useful  information.  An  example 
application  for  determining  resonant  behaviors  is  in  archeological  digging.  If  the  resonant 
length  of  a  partially  buried  fossil  is  known,  time  can  be  saved  in  the  process  of  unearthing 
the  object.  The  ability  to  discriminate  between  hollow  and  solid  targets  can  also  provide 
a  useful  tool  in  the  drug  war.  A  car  door,  which  is  supposed  to  be  hollow,  may  have  a 
different  acoustic  response  if  it  is  packed  with  drugs.  This  thesis  describes  the  techniques 
used  to  examine  several  different  targets  and  the  results  that  validate  our  approach. 

The  remainder  of  this  chapter  serves  as  a  review  of  acoustic-driven  techniques. 
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Chapter  2  presents  a  signal  processing  block  diagram  and  discusses  the  techniques  used  to 
analyze  the  data.  It  also  includes  an  air  propagation  model  proposed  to  improve  acoustic 
imaging  techniques  in  the  future.  Chapter  3  describes  the  experimental  setup  and  discusses 
the  necessity  to  build  an  anechoic  chamber  in  support  of  this  research.  Finally,  Chapter 
4  reveals  the  results  of  this  research,  and  concludes  with  what  future  work  may  include  in 
this  field. 


1.3  Review  of  Acoustic  Principles 

Before  delving  into  nonlinear  acoustic  detection,  a  review  of  the  fundamentals  of 
acoustics  is  necessary.  The  remainder  of  this  chapter  discusses  fundamental  wave  propaga¬ 
tion,  basic  acoustic  interactions,  and  what  assumptions  are  made  for  the  work  presented 
in  this  thesis.  This  will  provide  the  reader  with  appropriate  background  knowledge  to 
understand  the  physical  phenomena  that  take  place. 


1.3.1  Ideal  Linear  Wave  Equation 


The  ideal  linear  wave  equation,  which  is  presented  here,  is  derived  from  first  prin¬ 
ciples  by  Blackstock[3].  In  this  derivation,  the  concept  of  a  fluid  particle  is  introduced  in 
order  to  neglect  the  random  path  and  velocities  of  individual  molecules.  A  fluid  particle’s 
size  can  be  defined  as  the  smallest  volume  spanned  by  a  length  of  the  same  order  magnitude 
as  the  average  mean  free  path  of  a  molecule.  For  air,  this  length  l  is  approximately  0.1  mm 
and  the  number  of  molecules  contained  in  l3  is  about  25,000.  This  volume,  shown  in  Figure 
1.3  is  used  to  derive  the  equation  of  continuity.  In  Figure  1.3,  p  defines  the  density  of  the 
fluid,  u  defines  the  fluid  velocity,  and  S  defines  the  surface  area  of  the  flow.  The  net  change 
in  the  mass  of  this  fixed  volume  can  be  written  mathematically  as 

d 

—  ( SpAx )  =  puS \x  -  puS\x+Ax •  (1-1) 


Rearranging  terms  and  taking  the  limit  as  Ax  — ►  0  yields  the  equation  of  continuity, 


dp  d(pu) 
dt  dx 


(1.2) 


The  next  physical  property  to  model  is  the  conservation  of  momentum.  Figure  1.4 
models  the  momentum  flow  and  pressure  forces  acting  on  the  fluid  particle.  The  momentum 
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Figure  1.3:  Mass  flow  through  the  surface  of  a  cylindrical  fluid  particle[3]. 


Figure  1.4:  Momentum  inflow  and  outflow  and  forces  acting  on  cylindrical  fluid  particle[3]. 


inflow  and  outflow  as  well  as  the  pressure  acting  on  the  fluid  particle  can  be  modeled  using 


puSAx )  =  pu2S\x  -  pu2S\x+Ax  +  PS\X  -  PS\x+Ax, 


(1.3) 


where  P  is  the  pressure  acting  on  the  surface  of  the  fluid  particle.  Equation  (1.3)  can  be 
simplified  by  dividing  through  by  S' Ax  and  applying  Equation  (1.2).  Taking  the  limit  as 
Ax  — >  0  yields  the  conservation  of  momentum  equation, 


'  du  du\  dP 

p  |  —  +  u—  )  +  —  =  0. 


dt  '  “ dx  J  ’  dx  ^  ^ 

Finally,  the  ideal  linear  wave  equation  is  based  upon  a  thermodynamic  equation 
of  state  such  as  the  ideal  gas  law,  PV  =  nRT  ,  where  V  is  a  gas  volume,  n  is  the  number 
of  moles,  R  is  the  ideal  gas  constant,  and  T  the  absolute  temperature.  For  gases,  the  most 
frequently  used  isentropic  equation  of  state  is  the  adiabatic  gas  law, 


-W- 

PoJ  \Po 


(1.5) 
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where  7  is  the  ratio  of  specific  heats  and  po  and  po  are  the  static  values  of  P  and  p  respec¬ 
tively.  A  more  useful  equation  of  state  valid  for  both  liquids  and  gases  is  a  Taylor  series 
expansion  of  the  general  isentropic  equation  of  state  about  the  condensation  (p  —  po)/po  so 
that 


P  =  po  +  A 


P~  Po 
Po 


B 

+  2! 


P~  Po 
Po 


2+C(p-p0'3 


3! 


Po 


+ 


(1.6) 


The  coefficients  A,  B,  C...are  determined  from  experimental  observations  or  calculated 
using  nonlinear  approximations.  Equation  (1.6)  will  be  later  used  to  provide  a  third-order 
nonlinear  wave  equation  approximation.  Introducing  the  variable  for  sound  speed 


c2  = 


dP 


dp  ’ 


and  rearranging  terms,  Equation  (1.6)  becomes 


P  =  <W 


B  5p  C 

1  H - -  H - 

2\A  po  3!A 


A,  +... 

Po 


(1.7) 


(1.8) 


In  Equation  (1.8),  p  represents  the  excess  pressure  defined  by  p  =  P  —  po  and  5p  is  the 
excess  density  defined  by  bp  =  p  —  po.  Equations  (1.2),  (1.4),  and  (1.8)  may  be  made  linear 
by  using  a  small  signal  approximation.  This  approximation  is  valid  for  even  the  loudest 
sounds  that  are  experienced  on  a  day-to-day  basis.  The  small  signal  approximation  is  based 
on  the  assumption  that  excitation  pressures,  which  affect  density  and  particle  velocity,  are 
very  small  in  comparison  to  steady-state  values.  This  can  be  seen  by  the  fact  that  134 
dB  sound  pressure  level  (SPL),  which  is  equivalent  to  an  excitation  pressure  of  100  Pa,  is 
miniscule  in  comparison  to  101.325  kPa  ambient  air  pressure.  Using  the  assumption  that 
\p\  <C  Po£qi  we  can  linearize  Equations  (1.2),  (1.4),  and  (1.8)  to  yield 

(1.9) 


85  p  du 

~m+pWx  =  0' 


du  dp 
pdi  +  fric=  ’ 
P  =  C205p. 


(1.10) 

(1.11) 


Combining  Equations  (1.9-1.11)  produces  a  single  second  order  differential  equation  of 
pressure  with  respect  to  distance  and  time, 


2  d^p  d2p  _ 
C°  dx2  dt2 


(1.12) 


Equation  (1.12)  is  the  linear,  non-dissipative,  acoustic,  planar  wave  equation  and  has  a  real 
solution  of  the  form 


p(x,  t )  =  Pq  cos(c vt  —  kx) 


(1.13) 
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where  Po  is  the  excitation  amplitude,  u ;  is  the  frequency  in  radians,  and  k  =  oj/c  is  the  wave 
number.  This  solution  to  the  wave  equation,  as  well  as  those  presented  in  Sections  1.3.2 
and  1.3.3  are  for  one-dinrensional  plane  wave  propagation,  which  assumes  that  wave  fronts 
are  uniform  in  the  y-z  plane  and  located  in  the  acoustic  far  field  defined  by  [9] 

d 2 

(1.14) 

where  d  is  the  diameter  of  an  acoustic  source  and  A  =  c//  is  the  wavelength  of  the  acoustic 
signal. 


1.3.2  Lossy  Linear  Wave  Equation 


Dissipation  of  an  acoustic  signal  is  primarily  due  to  thermoviscous  losses  associated 
with  the  media  in  which  it  travels.  Dissipation  is  often  ignored  for  audible  frequencies  in 
air  because  the  acoustic  attenuation  coefficient,  a,  is  proportional  to  the  square  of  the 
frequency.  However,  as  frequency  increases  to  ultrasonic  levels,  the  dissipation  reaches 
levels  high  enough  to  limit  propagation  to  only  a  few  meters.  In  order  to  develop  a  wave 
equation  to  model  this  effect,  we  must  consider  the  equations  of  continuity,  momentum, 
and  state  [3].  Viscous  effects  do  not  alter  the  continuity  equation,  and  are  ignored  in  the 
linearized  equation  of  state  because  they  appear  as  higher  order  terms.  Only  the  equation 
of  momentum  is  changed  in  this  derivation,  adding  an  additional  term  to  the  right-hand 
side  yielding  [3] 


du  dp  • ~d2u 

/9°  dt  dx  ^  dx2  ’ 


(1.15) 


where  p  is  the  shear  viscosity  coefficient  and  V  =  4/3  +  pb/p-  Equations  (1.9),  (1.11),  and 
(1.15)  are  now  combined  to  yield  the  linear  lossy  wave  equation  of  pressure  with  respect  to 
distance  and  time, 


vV  d2pdp  d2p  1  d2p 
Cg  dx2  dt  dx2  Cq  dt 2 

where  v  =  p/po  is  the  kinematic  viscosity  coefficient.  Equation  (1.16)  has  a  real  solution  of 
the  form 


p(x ,  t)  =  P^e  ax  cos  (iot  —  kx ), 


(1.17) 


where  a  is  the  attenuation  coefficient  in  Nepers  per  meter.  The  attenuation  coefficient, 
a,  can  be  calculated  for  both  viscous  and  thermal  dissipation.  Typically,  attenuation  co¬ 
efficients  are  computed  for  viscous  effects  and  altered  heuristically  to  account  for  thermal 


effects.  Both  viscous  and  thermal  attenuation  are  represented  here  using 


OLv  — 

^ th 


VVLO2 

^T 

(7  —  1)huj2 
^PoCqCp 


(1.18) 

(1.19) 


where  k  is  the  heat  conduction  coefficient  and  Cp  is  the  specific  heat  at  constant  pressure. 
These  two  absorption  coefficients  can  be  summed  to  produce  the  total  absorption  coefficient 

7-1' 


&tv  — 


lo2v 
H  L 


4  PB_ 

3^  p 


+ 


Pr 


(1.20) 


where  Pr  is  the  Prandtl  number  and  Pb/p  =  0.61.  Equation  (1.20)  is  very  close  to  the 
classical  definition  of  acoustic  attenuation  given  by  [10]  [11] 


^ classical 


c o2v 


4  7  —  1 

-  +  1 - 

3  Pr 


(1.21) 


For  air  at  20°C,  Pr  =  0.711,  7  =  1.402,  v  =  15.11  x  10  6  m2 /s,  and  c  =  343  nr/s,  Equation 
(1.20)  yields  a  thernroviscous  absorption  coefficient  of  1.84  x  10-11  f2,  where  /  is  in  hertz. 


1.3.3  Lossy  Nonlinear  Wave  Equation 


When  amplitude  and  frequency  become  very  large,  small  signal  models  tend  to 
break  down  and  lose  accuracy  due  to  the  nonlinearity  of  air.  In  this  case,  a  nonlinear 
propagation  model  must  be  derived  to  account  for  the  formation  of  higher  order  harmonics. 
The  nonlinear  derivation  here  will  focus  on  the  state  equation  presented  in  Equation  (1.6), 
while  keeping  terms  of  order  greater  than  one.  Referring  back  to  the  original  state  equation 
given  in  Equation  (1.5),  and  taking  the  derivative  indicated  in  Equation  (1.7),  an  expression 
for  sound  speed  is  given  by 


c 


2 


7 P 
P 


7  RT. 


(1.22) 


It  is  important  to  observe  that  P,  p,  and  T  in  Equation  (1.22)  are  total,  non-static,  values, 
which  now  make  sound  speed  non-constant.  Combining  Equations  (1.5)  and  (1.22),  we  can 
express  P  and  p  in  terms  of  c  as 


27 

(7-1) 


P  =  Po 


(1.23) 
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Using  Equation  (1.23),  we  can  eliminate  P  and  p  from  the  continuity  and  momentum 
equations  and  rewrite  them  as 


Dc  7  —  1 

Dt+  2  CV'U  =  ° 

(1.24) 

Du  2 

(1.25) 

j-  + - -cV  •  c  =  0, 

Dt  7  —  1 

respectively.  For  plane  waves,  Equations  (1.24)  and  (1.25)  reduce  to 


dc  dc  7  —  1  du 

dt+Udx+  2  C9x“° 

(1.26) 

du  du  2  dc 

(1.27) 

dt  +Udx  +  7-lCdx~°- 

Equations  (1.26)  and  (1.27)  represent  a  second  order  system  that  can  be  reduced  by  only 
taking  into  consideration  the  forward  traveling  wave  component  represented  by 

g  +  to  +  wg-O,  (1.28) 


with  nonlinear  parameter  (3  =  (7+1) /2.  After  expanding  u  as  a  Taylor  series  and  performing 
extensive  algebraic  manipulations,  a  second  order  solution  to  Equation  (1.28)  of  the  pressure 
p  with  respect  to  x  and  t  can  be  written  as 


p(x,t)  =  Pq  (  sin (ut  —  kx)  +  ~/3Mkx  sin  2(cut  —  kx) 


(1.29) 


where  M  =  uq/cq.  It  is  clear  from  Equation  (1.29)  that  the  second  harmonic  amplitude  is 
directly  proportional  to  frequency,  amplitude,  and  propagation  distance  and  can  be  ignored 
if  all  three  are  small  enough.  If  the  Mach  number  is  expressed  with  respect  to  pressure, 
M  =  P0/(Z0c0)  ,  the  second  order  nonlinear  equation  may  be  written  as 


p(x,  t)  =  Pq  sin(u;f  —  kx)  + 


Pq(3ujx 


sin  2  (ut  —  kx), 


(1.30) 


where  Zq  is  the  characteristic  impedance  of  air.  Since  the  second  harmonic  is  dependent 
on  the  square  of  the  amplitude,  large  excitation  pressures  are  the  predominant  cause  of 
harmonic  distortion.  This  is  evident  at  a  rock  concert  where  music  sounds  distorted  due  to 
large  propagation  distances  and  extreme  amplification.  Combining  Equation  (1.30)  with  the 
results  obtained  in  Section  1.3.2  results  in  the  lossy  nonlinear  second  order  wave  equation 


p(x,t) 


P0e~aix  sin (ut  -  kx)  + 


P^Pujx_a2X  s{n2^t_kxy 


(1.31) 
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Higher  order  solutions  to  the  nonlinear  wave  equation  exist,  but  very  little  is  know  about 
higher  order  nonlinear  parameters  for  air.  If  air  is  treated  as  an  ideal  gas,  the  second  order 
nonlinear  parameter  is  /3  =  1.2.  In  Chapter  2,  a  third  order  nonlinear  air  model  is  derived 
using  a  perturbation  method. 


1.3.4  Modulation  of  Sound  by  Sound 


This  section  discusses  the  physical  mechanisms  involved  in  using  ultrasonic  para¬ 
metric  arrays  for  standoff  analysis  of  targets.  In  1982,  Yoneyama[4]  discussed  the  nonlinear 
interaction  of  ultrasound  with  air  as  the  “scattering  of  sound  by  sound”.  This  work  is 
primarily  based  on  the  nonlinear  wave  equation  derived  by  Westervelt[12]  to  account  for 
deficiencies  of  classical  models  at  high  frequencies.  Westervelt’s  answer  to  the  high  fre¬ 
quency  problem  is  the  addition  of  sum  and  difference  frequencies,  which  accounts  for  the 
nonlinear  interaction  between  air  molecules  and  is  expressed  as 


w2  1  d  Ps 

V  Ps  ~2  5f2 


=  Po 


~'0 


dq 

dt: 


q  = 


p  d 


P  oco 


Pi- 


(1.32) 


In  Equation  (1.32),  ps  is  the  secondary  wave  pressure,  pi  is  the  primary  wave  pressure,  and 
P  is  the  second  order  nonlinear  parameter.  The  nonlinear  interaction  of  air  molecules  is  due 
to  the  difference  in  the  forces  in  the  compression  and  tension  phases  of  acoustic  propagation. 
At  any  point  along  its  propagation  axis,  a  longitudinal  acoustic  wave  produces  areas  of  high 
and  low  pressure.  In  the  high-pressure  (compression)  phase  of  propagation,  two  neighboring 
air  molecules  move  together  with  minimal  compression.  In  the  low-pressure  (tension)  phase 
of  propagation,  the  retreating  air  molecule  has  a  lesser  pull  on  its  neighbor  due  to  the 
weak  attraction  between  them.  There  is  slight  delay  before  the  molecule  being  pulled  on 
can  catch  up,  producing  a  “slapping”  or  cavitation  effect.  This  produces  a  rectification  at 
higher  frequencies  just  as  a  diode  does  in  an  electrical  circuit. 

This  natural  rectification  allows  the  demodulation  of  an  AM  waveform.  In  the  sim¬ 
plest  interactions,  two  pure  tones  f\  and  fi  mix  to  produce  second  order  inter-modulation 
products  of  f\  ±  f‘2  ■  Yoneyama  elaborated  on  this  simple  phenomenon  to  demodulate  an 
AM  waveform  with  secondary  pressure  given  by 


1  +  mg 


e~ax  sinw 


Ps(x,t)  =  pi 


(1.33) 
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where  m  is  the  modulation  index  and  g  represents  an  arbitrary  base  band  signal.  As  in 
Equation  (1.32),  ps  is  the  secondary  pressure,  which  in  this  application  seems  to  emanate 
from  thin  air.  Furthermore,  because  the  modulating  waveform  g  is  carried  by  a  high  fre¬ 
quency  signal  with  wavelength  much  smaller  than  the  size  of  the  transducer,  these  signals 
are  highly  directional  as  shown  in  Figure  1.5.  The  only  significant  drawback  of  using  the 


Figure  1.5:  Directivity  of  secondary  wave  at  10.0  kHz,  for  a  point  at  4  m,  m  =  0.5,  and 
input  voltage  of  10  V  [4], 


nonlinear  effect  of  air  to  demodulate  ultrasonic  signals  is  the  poor  efficiency  of  the  demod¬ 
ulation  process.  Another  way  of  expressing  Equation  (1.33)  is 

,3P2A  d2 


Ps(t)  = 


IQi: pqCqXol  dt2 


E 


x 

t - 

co 


(1.34) 


Here  A  is  the  cross  sectional  area  of  the  transducer  and  E(t)  is  the  modulation  envelope. 
Inspecting  Equation  (1.34),  it  can  be  seen  that  the  secondary  pressure  is  proportional  to 
the  second  derivative  of  the  square  of  the  modulation  envelope.  The  second  derivative  term 
produces  a  slope  in  the  frequency  domain  of  12  dB  per  octave,  and  the  square  term  adds 
significant  distortion  in  double-sideband  AM  modulation.  As  an  example,  using  Equation 
(1.34)  with  / 3  =  1.2,  A  =  0.2  m2,  a  =  0.7,  cq  =  348  m/s,  and  po  =  1.18  kg/m3,  a  130  dB 
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carrier  modulated  with  a  1  kHz  signal  produces  about  66  dB  of  audible  sound  at  1  nr.  With 
the  ability  to  only  place  sound  intensities  in  the  60  dB  range  on  the  surface  of  the  target,  it 
is  apparent  why  high  resolution  spectral  techniques,  like  those  discussed  in  Chapter  2,  are 
necessary  for  this  type  of  detection. 

1.3.5  Decibel  Scale 

The  time  is  taken  here  to  discuss  logarithmic  scales  as  they  are  used  in  the  field 
of  acoustics.  This  is  a  necessary  evolution,  as  it  helps  to  understand  the  relevant  intensities 
that  are  invoked  in  a  comparison  to  Ohm’s  law  as  it  relates  to  electricity.  Table  1.1  relates 
the  basic  measurable  acoustic  parameters  to  their  equivalent  electrical  properties. 


Table  1.1:  Conversion  between  acoustic  and  electric  variables 


Acoustic  Variable 

Electric  Variable 

Pressure  (P)  Pascals 

Voltage  (V)  Volts 

Particle  Velocity  (u)  m/s 

Current  (I)  Amps 

Acoustic  Impedance  (Z)  Pa  ■  s/m 3 

Resistance  (R)  Ohms 

Sound  Intensity  (I)  W/m2 

Power  (P)  Watts 

This  is  helpful  in  understanding  the  derivation  of  sound  pressure  level  (SPL)  as 
a  means  of  representing  acoustic  power.  Using  the  relationship  between  all  four  acoustics 
parameters  in  Table  1, 


7=p=L= H 

u  v?  I  : 


(1.35) 


it  is  easy  to  see  that  sound  intensity  is  the  same  as  electric  power.  Just  as  logarithms  are 
used  to  express  electric  power  with  respect  to  a  reference  power,  they  are  used  to  express 
sound  intensity  with  respect  to  a  reference  intensity  of  20  /r Pa .  This  is  considered  to  be 
the  lowest  intensity  perceivable  by  the  human  ear  at  2  kHz.  Representing  sound  intensity 
in  logarithmic  scale  can  be  accomplished  by  either  of  the  following  equations, 


SPLdB  =  101og10  =  201og10  ■  (1-36) 

This  representation  is  useful  in  measuring  the  wide  range  of  amplitudes  that  are  typical  in 
everyday  environments.  Table  1.2  presents  several  examples  of  sound  pressure  level  with 
respect  to  excitation  amplitudes  from  various  sources. 
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Table  1.2:  Examples  of  sound  pressure  and  corresponding  sound  pressure  level  (SPL) 


Source 

Sound  Pressure 

Sound  Pressure  Level 

1  atmosphere 

101325  Pa 

191.085  dB 

Thermoacoustic  device 

12000  Pa 

176  dB 

Jet  engine  at  30  nr 

630  Pa 

150  dB 

Threshold  of  pain 

100  Pa 

130  dB 

Jack  hammer  at  1  nr 

2  Pa 

100  dB 

Hearing  damage 

6  x  10"1  Pa 

85  dB 

Passenger  car  at  10  nr 

2  x  10~2  -  2  x  10_1  Pa 

60-80  dB 

Normal  talking  at  1  nr 

2  x  10~3  -  2  x  10~2  Pa 

40-60  dB 

Very  calm  room 

2  x  10~4  -  6  x  10“4  Pa 

20-30  dB 

Auditory  threshold  at  2  kHz 

20  x  10~e  Pa 

0  dB 

1.3.6  Acoustic  Absorption  and  Reflection 

Because  this  thesis  addresses  the  interaction  of  ultrasonic  acoustics  on  solid  targets, 
a  discussion  of  the  interaction  of  sound  at  a  boundary  is  included.  Just  as  in  transmission 
line  theory,  when  a  voltage  traveling  on  a  line  with  characteristic  impedance  Zq  encounters 
a  circuit  with  a  different  characteristic  impedance  Zi,  some  energy  is  reflected  while  the 
load  absorbs  the  rest.  Ultrasonic  techniques  for  non-destructive  testing  were  first  employed 
by  Sokoloff[13]  in  1929,  and  have  become  commonplace  today  with  the  growth  of  microelec¬ 
tronics  and  high  sensitivity  transducers.  However,  due  to  the  large  impedance  mismatch 
between  air  and  solids,  most  ultrasonic  inspection  techniques  require  the  test  item  to  be  im¬ 
mersed  in  water  or  use  a  gel  coupling  medium[5].  Neither  of  these  methods  are  suitable  for 
standoff  detection  and  until  recently,  technology  did  not  exist  to  efficiently  transmit  ultra¬ 
sound  for  non-contact  measurements[14].  With  the  recent  advent  of  piezo-electric  ceramic 
transducers  incorporating  a  polymer  matching  layer,  parametric  arrays  can  be  constructed 
to  provide  highly  directional,  high  amplitude  ultrasound. 

In  this  derivation,  only  plane  waves  of  normal  incidence  will  be  considered  to 
avoid  rigorous  trigonometric  equations.  The  incident,  reflected,  and  transmitted  values  for 
pressure,  intensity,  and  particle  velocity  will  be  defined  by  subscripts  i,  r,  and  t  respectively 
and  can  been  seen  in  Figure  1.6.  Just  as  in  electromagnetic  theory,  the  characteristic 
impedance  describes  the  conductive  properties  of  a  medium  and  is  defined  for  air  by 

Zq  =  PqCq.  (1-37) 


Conservation  of  energy  requires  the  use  of  boundary  conditions  at  x  =  0  to  define  the 
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incident  w 

Pi  1  i  V, 

x=0 

'ave  medium  1 

=  Pici 

reflected  wave 
'  Pr  lr  Vr 

acoustic  boundary 

X 

transmitted  wave 
’  P,  1,  v, 

medium  2 

^2  =  P2P2 

Figure  1.6:  Transmission  and  reflection  of  an  acoustic  wave  at  normal  incidence  at  a  bound¬ 
ary  of  two  materials  with  different  impedance  [5] . 

forward  and  backward  traveling  wave  with  pressure  pr,  particle  velocity  ur,  and  intensity 
Ir  defined  by 


U%  T  Ur  —  Ut\x=01 

Pi+Pr  =  Pt\x=0,  (1.38) 

b  Ir  —  lt\x=0- 

Solving  for  the  incident,  reflected,  and  transmitted  pressures  with  respect  to  density,  sound 
speed,  and  particle  velocity  yields 


Pi  —  Pi  Cl  U i  |  x=0  > 

Pr  —  PlClUr\x=ih  (1.39) 

Pt  =  P2C2Ut\x=0- 

By  combining  Equations  (1.39)  and  (1.37),  Equation  (1.38)  becomes 


Ui  +  Ur  —  Ut  |  oi=0  5 
UiZ\  UrZ\  —  'U’fZ2\x=0 


(1.40) 
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and 


Pi  Pr  —  Pt\x=Oi 
Pi_  _  Pr_  _  Pt_\ 

ry  ry  ry  \  X — 0* 

Zj\  Zj  i  Zj2 


(1.41) 


If  we  define  the  reflection  coefficient  as  the  ratio  of  the  reflected  wave  to  the  incident  wave. 
Equations  (1.38),  (1.40),  and  (1.41)  can  be  combined  to  get 


R 


p  ~ 


Pr 

Pi 


Z'2  ~  Z\ 
Z-2  +  Z\ 


(1.42) 


The  transmission  coefficient  is  defined  in  a  similar  fashion  and  is  computed  using 


Pt  _  2Zi 
Pi  Z‘2  +  Z\ 


(1.43) 


Table  1.3  lists  the  acoustic  impedance  of  several  common  materials  and  the  reflection  and 
transmission  coefficient  for  each  one  with  respect  to  air.  The  values  listed  help  further 
clarify  why  a  high-resolution  spectral  technique  is  needed  to  analyze  resonant  responses  in 
this  standoff  nonlinear  acoustic  detection  technique. 

Table  1.3:  Reflection  and  transmission  coefficients  for  some  common  materials  with  respect 
to  air,  Zair  =  413.5  MRayl. 


Material 

Acoustic 

Impedance  (Rayl) 

Reflection 

Coefficient 

Transmission 

Coefficient 

Air 

413.5 

0 

1 

Water 

1.48  x  106 

0.99944 

0.00055 

Wood  (pine) 

1.57  x  106 

0.99947 

0.00052 

Fiberglass 

2.86  x  106 

0.99971 

0.00028 

Concrete 

8.00  x  106 

0.99989 

0.00010 

Glass 

12.5  x  106 

0.99993 

0.00006 

Sand 

19.7  x  106 

0.99995 

0.00004 

Steel 

46.0  x  106 

0.99998 

0.00001 
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Chapter  2 

Signal  Processing 

2.1  Introduction 

This  chapter  focuses  on  the  development  of  signal  models  and  signal  processing 
algorithms  to  overcome  the  low  signal  to  noise  ratio  (SNR)  discussed  in  Chapter  1.  Due 
to  the  high  loss  of  the  nonlinear  demodulation  of  ultrasonic  signals  and  the  poor  acoustic 
coupling  to  the  target,  reflected  signals  are  very  close  to  the  noise  floor.  For  the  anechoic 
chamber  used  in  these  experiments,  the  noise  floor  associated  with  acoustic  noise  plus 
measurement  noise  is  approximately  15  dB.  However,  the  noise  floor  and  SNR  are  closely 
related  to  the  sampling  frequency  and  sample  size.  If  measurements  are  made  using  a 
sampling  rate  of  500  kHz  and  taking  500,000  samples,  a  noise  floor  of  6-8  dB  is  achievable, 
but  highly  impractical  due  to  computation  time. 

What  is  required  is  a  high-resolution  technique  that  is  capable  of  analyzing  rela¬ 
tively  small  data  samples.  The  goal  is  to  sweep  the  difference  frequency  in  sufficiently  small 
increments  to  observe  narrowband  phenomenon  in  the  frequency  response.  Unfortunately, 
the  frequency  resolution  is  inversely  proportional  to  the  SNR  due  to  the  large  size  of  the 
data  set.  For  example,  if  we  want  a  frequency  sweep  from  0-10  kHz  in  100  Hz  increments, 
and  have  a  sample  rate  of  200  kHz  with  N  =  20,  000  samples,  our  data  set  has  2  million 
measurements.  It  is  easy  to  see  how  quickly  data  sets  can  grow  if  more  frequency  resolution, 
bandwidth,  or  SNR  is  needed. 


17 


2.2  Third  Order  Nonlinear  Air  Model 

Before  signal  processing  can  be  performed  on  any  data  set,  a  comprehensive  data 
model  must  be  established  to  define  the  signal  to  be  analyzed.  For  example,  the  data  model 
using  Fourier  analysis  defines  the  data  to  be  periodic  in  time.  Data  models  also  help  predict 
how  a  signal  varies  to  changes  in  the  medium  which  it  propagates  in.  An  example  directly 
related  to  this  research  would  be  the  change  in  sound  speed  with  respect  to  temperature 
and  pressure,  which  is  the  cause  of  its  nonlinear  behavior.  Scientists  have  used  nonlinear 
analysis  to  characterize  microwave  devices,  biological  tissues,  organic  liquids,  and  especially 
metal  components  [15,  16,  17,  18].  Nonlinear  characterization  has  been  especially  useful 
in  characterizing  small  bubbles  and  defects  within  living  tissue  and  metal  castings.  The 
underlying  theory  points  to  the  fact  that  defects  and  air  pockets  significantly  impact  second 
order  nonlinear  harmonic  generation.  Thus  a  measurement  of  the  nonlinearity  of  a  material 
compared  to  what  it  is  known  to  be  can  detect  defects  that  are  not  visible  using  other 
methods.  Chen[16]  reports  that  higher  harmonics  (n  >  3)  possess  superior  characteristics 
for  improved  clarity  and  sharper  contrast  in  biological  imaging. 

Third  order  nonlinear  models  have  been  developed  for  biological  tissues,  metals, 
and  microwave  devices  for  both  their  academic  understanding  and  prospective  applications 
[16,  19].  Nonlinear  behavior  of  acoustics  in  air,  especially  at  ultrasonic  frequencies,  have 
been  of  interest  to  scientists  for  quite  some  time.  However,  due  to  the  lack  of  advanced 
measurement  equipment  such  as  piezoelectric  microphones  and  transducers,  nonlinear  ul¬ 
trasonic  acoustics  has  mostly  remained  a  theoretical  interest.  With  the  advent  of  matching 
layers  to  transmit /capture  a  significant  amount  of  energy  from  a  transducer,  the  field  of 
ultrasonic  acoustics  has  experienced  a  resurgence  of  activity.  Both  the  fields  of  biology  and 
microwave  technology  understand  the  value  of  accurate  behavioral  models,  which  should 
carry  over  to  the  field  of  acoustics.  Furthermore,  the  nonlinear  model  presented  in  this 
section  goes  beyond  the  treatment  of  air  as  an  ideal  gas,  and  the  limitations  imposed  by 
other  third  order  models. 

2.2.1  Bessel  Approximation 

A  classical  treatment  of  nonlinear  acoustics  is  the  Bessel-Fubini  solution,  which 
allows  one  to  trace  the  growth  of  higher  order  harmonics.  The  modeling  of  acoustic  waves 
may  be  done  with  the  parameters  of  pressure  ( p ),  particle  velocity  (u),  and  displacement 
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(£).  In  this  treatment,  we  will  consider  the  parameter  of  particle  velocity,  and  represent  the 
ideal  linear  wave  equation,  Equation  (1.13)  with  respect  to  particle  velocity  as[20] 

u(x,  t)  =  uq  cos(ut  —  kx).  (2-1) 


If  we  define  sound  speed  with  respect  to  particle  velocity  as  [5] 

7  +  1 

C  =  C0  H - —  u, 


(2.2) 


and  consider  the  second  order  parameter  of  nonlinearity  f3,  we  can  describe  the  distortion 
of  a  sinusoidal  waveform  mathematically  as 

x 


u(x,  t)  =  u o  sin  uj  i  t  — 


0  <  x  <  X£>. 


(2.3) 


co  +  /3u/ 

In  Equation  (2.3),  xd  represents  the  discontinuity  distance,  or  the  distance  at  which  the 
waveform  changes  from  a  sinusoidal  nature  to  an  almost  sawtooth  shape.  This  phenomenon 
is  caused  by  the  difference  in  wave  propagation  speed  between  the  crests  and  troughs.  The 
crests  represent  high  pressure  (high  density)  areas  while  the  troughs  represent  low  pressure 
(low  density)  areas.  Since  sound  travels  faster  in  higher  density  media,  the  crest  of  the 
wave  tries  to  overtake  the  trough,  forming  a  sawtooth  shape  as  can  be  seen  in  Figure  2.1. 
The  distance  at  which  this  occurs  is  based  on  the  frequency  and  amplitude  of  the  excitation 
waveform  and  is  represented  as 

Xd  ~  u  a  *  r  ’  (2-4) 


k/3M 


where  M  is  the  acoustic  Mach  number  M  =  uq/cq.  The  frequency  content  of  Equation 
(2.3)  can  better  be  seen  by  representing  it  as  a  Fourier  series: 


u(x,  t)  =  uq  Bn  sin (n(u+  —  kx)),  x  <  xd 


Bn  =  - 


1 

7 T  -'0 


71=1 
2tt 


/*Z7T 

/  sin(cuf  —  kx  +  (3)  sin (n(ixt  —  kx))d{ut  —  kx). 

Jo 


The  integration  of  Equation  (2.6)  yields 


Bn  =  2^Jn  (—)  , 
nx  \xd  J 


(2.5) 

(2.6) 

(2.7) 


for  nth  order  Bessel  functions,  Jn,  of  the  first  kind.  Putting  Equation  (2.7)  into  Equation 
(2.5),  one  obtains  a  final  solution  with  respect  to  pressure  of [5,  20] 


p(x,  t)  =  2 P0 


oo  jim. 
—  Jn\xD 


n= 1 


nx 
X  D 


sin  (n(iot  —  kx)), 


x  <  xd- 


(2.8) 
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Figure  2.1:  Image  of  a  30  kHz,  117  Pa  sine  wave  being  distorted  by  nonlinear  effects  of  the 
air.  Solid  line  represent  propagation  distance  x  =  0,  dotted  line  x  =  4  nr,  and  dashed  line 
x  =  8  nr  is  the  theoretical  waveform  shape  beyond  discontinuity  distance  x  =  xe>. 


This  equation  describes  the  spectral  content  of  the  higher  order  harmonics  as  their  ampli¬ 
tudes  grow  at  the  expense  of  the  fundamental  frequency.  Notice  however  that  this  model 
only  works  for  distances  less  than  the  discontinuity  distance,  xr>.  This  severely  limits  the 
usefulness  of  this  approximation  for  high-frequency,  high-amplitude  signals  since  the  dis¬ 
continuity  distance  may  be  very  small. 

2.2.2  Derivation  of  Nonlinear  Parameters 

One  solution  to  providing  a  more  accurate  model  at  further  ranges  than  the  Bessel- 
Fubini  equation  is  called  the  perturbation  method.  This  approach  begins  with  the  nonlinear 
equation  of  state  repeated  here  for  convenience,  [3] 

(^)+§(^)2+i(— 


P  =  Po  +  A 


Po 


(2.9) 
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In  Equation  (2.9),  B/A  and  C/A  are  the  ratio  of  the  quadratic  to  linear  terms  and  cubic 
to  linear  terms  of  the  Taylor  series  expansion.  Their  definitions  can  be  represented  by 
partial  derivatives  of  pressure  with  respect  to  density,  evaluated  at  equilibrium  conditions 
and  constant  entropy,  which  are  denoted  by  subscripts  “0”  and  “s”  [16] .  The  first,  second, 
and  third  order  terms  of  Equation  (2.9)  can  be  solved  for  using 


A  =  po 
B  =  pl 

c  =  pI 


dp\ 

dp) 

s,0 

(2.10) 

d2p 

dp2 

)  s,0 

(2.11) 

d3p' 

dp3 

)  s,0 

(2.12) 

Noticing  that  (dp/ dp)  =  Cq  in  Equation  (2.10),  the  equation  for  A  can  be  simplified.  Taking 
the  ratio  of  B/A  and  C/A  further  reduces  the  complexity  to  only  second  order  derivatives 
of  sound  speed  with  respect  to  pressure  given  by [21] 


A  =  p0cl 
B 

—  -  2p0c0 

C_3(B 
A  ~  2  \A 


dc\ 
dp)  s,  o 

+  2  p2c3 


(2.13) 

(2.14) 

(2.15) 


The  problem  with  evaluating  Equations  (2.14)  and  (2.15)  is  finding  a  nonlinear  equation 
for  sound  speed  in  air.  Typically,  air  is  taken  as  an  ideal  gas,  and  sound  speed  is  linearized 
to  \J %BT JM.  This  problem  may  be  overcome  by  rearranging  Equation  (1.23),  which  is 
derived  from  the  adiabatic  state  equation  and  definition  of  sound  speed,  to  be 

c  =  c0(-)2\  (2.16) 

\PoJ 

This  definition  for  sound  speed  is  infinitely  differentiable  and  easily  provides  a  solution  to 
Equations  (2.14)  and  (2.15).  The  first  and  second  derivatives  of  Equation  (2.16)  are  given 
by 


dc  1 

(£) 

1  (7  -  l)co 

dp 

27  P 

d2c 

(£)*h -i)>« 

(£) 27  (t,-1)co 

dp2 

4^2p2 

2  7P2 

(2.17) 

(2.18) 
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These  derivatives  produce  solutions  to  Equation  (2.9)  of  B/A  =  .3998  and  C/A  =  —.2398, 
which  are  very  close  to  the  values  for  an  ideal  gas  of  B/A  =  ”/  —  1  and  C/A  =  7  —  2,  where 
7  =  1.4  for  air.  It  is  important  to  notice  from  Equations  (2.17)  and  (2.18)  that  the  nonlinear 
parameters  (3  =  1  +  \  ^  and  rj  =  1  +  5  §  are  proportional  to  the  excitation  pressure,  which 
is  in  keeping  with  the  results  of  Chen  and  Keller. 


2.2.3  Perturbation  Analysis 


Using  the  results  for  the  nonlinear  parameters  /3  and  ? 7  obtained  from  Section  2.2.2, 
we  can  apply  a  perturbation  solution  to  solve  for  the  amplitudes  of  higher  order  harmonics. 
Following  the  derivation  given  by  Chen,  we  can  model  a  finite  amplitude  plane  wave  in  a 
lossless  medium  using  [16] 


dp  dp  du 

m+U&c+Pfrc  =  ° 
du  du  dp 

Pm+pUd^+d-X=° 

P  =  P{P,  s )• 


(2.19) 


It  is  possible  to  find  a  solution  to  Equation  (2.19)  by  expressing  p,  p,  and  rasa  series, 

OO 

p-iy*1 

n= 0 
00 

u  =  J2u^n) 

n= 0 

OO 

p  =  YJp(n\  (2.20) 

n= 0 


and  substituting  these  back  into  Equation  (2.19).  We  can  now  independently  express  the 
amplitudes  of  the  second  and  third  harmonics  as  second  order  partial  differential  equations 
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of  pressure  with  respect  to  distance  and  time,  and  given  by 


d2p^  i  o2p (2) 

dx2  Cq  dt 2 


d_ 

dx 


(p0  lco(P  ~  !) 


d(p^)2 

dx 


du^ 

dt 


+  Pou{1) 


+Pq  W  —  1) 


d2(p^)2 

dx2 


+ 


d2p^u^ 

dxdt 


Po  lco (P  ~  !) 


df*(pW)2 

dt 2 


(2.21) 


<92p( 3) 

dx2 


hrsr  =  -2pP{li  ~  ~  ^  ~  1)fe(',(1>)3 


f  p<» 


+ 


d_ 
dx 

\ 

a2 

dtdx 


duM 

~dt 


+  P 


(2)du(1'>  ( X)du ^  n-,  (2)®u^ 

{2)  -w—  +  Pquw  —r —  +  p(1)u(1)  +  pqvS2)  - r— 

dt  dx  dx  dx 


(pWu<V  +  pWuW). 


(2.22) 


In  Equations  (2.21)  and  (2.22),  the  superscript  in  parenthesis  refers  to  the  order  of  the 
harmonic  and  is  not  an  exponent.  Assuming  the  wave  at  x  =  0  is  sinusoidal,  Equations 
(2.21)  and  (2.22)  can  be  solved  with  respect  to  the  pressures  p^  and  p ^  respectively. 
Summarizing  the  total  solution,  we  can  represent  the  third  order  nonlinear  propagation 
model  as 


p^\x,t)  =  Poe  “l3;sin(u;t  —  kx) 

(2.23) 

p(2\x,t)  =  e~a2X  sin(2((jf  —  kx)) 

2/Jocg 

(2.24) 

p^(x,t)  =  GfPo(Ff/32x2  +  2px)e~a3X  sm(3(ujt  —  kx)) 

(2.25) 

G=  25 

2P6co 

(2.26) 

67T 

J  1 

CO 

(2.27) 

where  the  e~ax  term  has  been  added  to  each  harmonic  to  account  for  dissipative  losses  in 
a  medium.  The  individual  harmonics  can  then  be  summed  to  produce  the  final  solution  of 

3 

p(x,t)  =  ^p(x,t){n).  (2.28) 

n=  1 


Figure  2.2  shows  how  the  waveform  changes  as  it  propagates  in  the  positive  x  direction. 
Notice  how  the  approximations  all  look  similar  until  the  discontinuity  distance  xjj,  where 
the  Bessel-Fubini  solution  breaks  down.  At  this  point,  the  perturbation  and  ideal  gas 
solution  become  more  accurate  as  the  amplitude  of  the  Bessel-Fubini  approximation  falls 
too  low.  These  plots  look  time  reversed  from  those  in  Figure  2.1  because  they  are  plotted 
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with  respect  to  time  and  not  distance,  thus  flipping  them  left  to  right.  Figure  2.3  shows  a 
more  detailed  view  of  how  the  amplitudes  of  the  first,  second,  and  third  harmonics  of  the 
perturbation,  Bessel-Fubini,  and  ideal  gas  approximations  grow  and  decay  with  distance. 
Finally,  Figure  2.4  shows  how  the  total  acoustic  power  decreases  with  distance.  These  plots 
can  be  reproduced  and  altered  by  using  the  “Air_Model_GUI.m”  program  given  in  Appendix 
C. 


Figure  2.2:  30  kHz  sine  wave  at  120  dB  SPL  amplitude  plotted  at  x  =  0,  x  =  xd,  and 
x  =  2xd-  Notice  how  the  Bessel-Fubini  approximation  falls  away  beyond  x  =  xd,  which  is 
where  this  model  becomes  invalid. 


2.2.4  Fourier  Transform  Analysis 

The  Fourier  Transform  was  initially  used  in  our  spectral  analysis,  but  was  quickly 
dismissed  due  to  the  necessity  of  long  data  sets  to  achieve  sufficient  SNRs.  However,  the  use 
of  the  Discrete  Fourier  Transform  (DFT)  was  an  integral  part  of  our  processing  algorithm, 
in  its  role  of  implementing  a  prewhitening  filter.  On  this  account,  we  provide  here  a  brief 
discussion  on  its  application. 
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1st  Harmonic  2nd  Harmonic  3rd  Harmonic 


Distance  (x)  meters 


Figure  2.3:  Harmonic  amplitudes  of  30  kHz  sine  wave  at  120  dB  SPL  amplitude.  Notice 
how  all  three  amplitudes  are  reasonably  close  for  x  <  xd  =  7.6  m,  which  is  where  the 
Bessel-Fubini  approximation  becomes  invalid  and  the  perturbation  approximation  becomes 
more  accurate. 


Time  domain  signals  of  incident  and  reflected  pressure  waveforms  are  fundamental 
to  a  pursued  detection  procedure.  While  the  incident  waveform  does  not  contain  any  data 
about  the  composition  of  the  target,  it  does  provide  a  reference  energy  level  as  a  means  of 
comparing  the  amount  of  acoustic  energy  absorbed  by  solid  and  hollow  targets  of  the  same 
shape.  The  incident  waveform  also  confirms  that  the  parametric  array  is  correctly  pointed 
at  the  target,  ensuring  that  it  is  receiving  maximum  power.  In  order  to  maintain  the  best 
spectral  resolution  from  a  DFT,  we  must  preserve  the  tinre-bandwidth  product [22] 

TeBe  =  1.  (2.29) 


In  Equation  (2.29), 


1  _  1  _  Fs 
Te  NT  N  ’ 


(2.30) 


where  Te  is  the  observation  interval,  Be  is  the  frequency  resolution,  N  is  the  number  of 
samples,  and  Fs  the  sampling  rate.  By  capturing  50,000  samples  at  a  sample  rate  of  50 
kHz,  we  are  able  to  achieve  a  frequency  resolution  of  1  Hz,  which  is  the  maximum  achievable 
by  the  signal  generators  being  used.  The  calculated  frequency  resolution  is  then  used  in 
generating  the  frequency  “bins”  in  the  single  sided  DFT  using  [22] 


X[k] 


N—l 

—  V  x[n}e~j2nk% 
N  '  L  J 

72—0 


(2.31) 
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Total  Acoustic  Intensity 


5  10 

Distance  (x)  meters 


Figure  2.4:  Total  acoustic  intensity  of  30  kHz  sine  wave  at  120  dB  SPL  amplitude  plotted 
with  respect  to  distance.  Notice  how  all  propagation  models  start  at  120  dB  at  x  =  0,  and 
how  the  Bessel-Fubini  model  becomes  more  inaccurate  beyond  x  =  xd- 

and  mapped  to  the  frequency  domain  using  the  substitution 

X[f]  =  X[k] \k=m  (2.32) 

Fs 

Equation  (2.31)  is  implemented  as  a  MATLAB  function  call  “generaLft.m”  and  a  “mexJft.nr” 
executable  to  improve  run  time.  The  code  for  these  functions  is  available  in  Appendix  A. 
Also  necessary  in  implementing  the  prewhitening  filter,  is  an  Inverse  Discrete  Fourier  Trans¬ 
form  (IDFT)  to  reconstruct  the  filtered  time  domain  signal.  This  IDFT  takes  the  real  and 
imaginary  Fourier  coefficients  and  produces  a  time  domain  representation  while  preserving 
the  phase  information.  If  the  phase  information  is  neglected  in  the  reconstruction  process, 
an  accurate  time  domain  signal  is  still  produced  with  respect  to  magnitude.  This  provides 
an  ability  to  build  “perfect”  filters  in  the  frequency  domain  and  transform  them  to  the  time 
domain  for  convolution  with  an  input  signal.  Reconstruction  of  a  time  domain  signal  from 
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the  single  sided  DFT  is  computed  using 

N- 1 

x[n]  =  ^  X[k\ej2*k% ,  (2.33) 

k= 0 

where  X[k\  is  a  complex  Fourier  coefficient.  Since  the  DFT  and  IDFT  are  implemented  in 
C,  all  values  must  be  “real”.  In  order  to  reconstruct  the  original  using  only  real  values, 
Equation  (2.33)  can  be  rewritten  as 

iV-l 

x[n]  =  ^  Re(X[k])  cos(2irk^)  —  Im(X[k])  sin(27r£:-^).  (2.34) 

k= o 

Equation  (2.34)  is  implemented  as  a  MATLAB  function  call  “generaLift.m”  and  a  “mexdfFt.m” 
executable  as  well.  Code  for  these  functions  is  available  in  Appendix  B.  Figure  2.1  (A)  shows 
an  input  signal  with  frequency  domain  representation  (B),  perfect  reconstruction  (C),  and 
no- phase  reconstruction  (D). 


(A)  Input  Signal  (2  kHz  Sinusoid) 


2 


(B)  DFT  using  "mex_fft.c" 


(C)  Reconstructed  with  "mexjfft.c" 


Time  (milliseconds) 


Figure  2.5:  Demonstration  of  “mex_fft.c”  and  “mexjfft.c”  with  2  kHz  input  signal  (A),  fre¬ 
quency  domain  representation  (B),  perfect  reconstruction  (C),  and  no-phase  reconstruction 
(D).  ‘ 


The  Multiple  Signal  Classification  (MUSIC)  algorithm  was  chosen  due  to  its  su¬ 
perior  performance  in  comparison  to  the  DFT.  However,  the  MUSIC  algorithm  makes  two 
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assumptions  about  the  input  data.  First,  it  assumes  that  the  data  is  zero  mean.  Second, 
it  assumes  that  the  additive  noise  is  white.  This  initially  presented  a  problem,  as  the  mea¬ 
surement  equipment  that  is  used  generates  1/f  “pink”  noise,  otherwise  known  as  flicker 
noise  when  referring  to  electronic  devices[23].  Figure  2.2  depicts  the  signal  processing  block 
diagram  that  was  developed  to  condition  the  data  prior  to  implementing  the  MUSIC  al¬ 
gorithm.  The  first  operation  that  is  performed  on  the  data  set  is  to  average  it  so  that 


Figure  2.6:  Signal  processing  block  diagram 

its  mean  is  zero.  Since  our  signals  are  theoretically  purely  sinusoidal,  this  step  is  largely 
precautionary.  However,  if  there  is  a  DC  component  added  by  the  measurement  equipment, 
this  step  will  eliminate  it.  Averaging  is  simply  and  mathematically  described  by 

x[n]avg  =  x[n\  -  E{x[n}},  (2.35) 

where  E  represents  the  expected  value  of  a  one  dimensional  vector  (time  series).  The  input 
signal  applied  to  the  prewhitening  filter,  thus  will  have  white  noise. 

2.2.5  Prewhitening  Filter 

As  previously  stated,  the  MUSIC  algorithm  assumes  that  the  time  domain  input 
signal  is  zero  mean  and  contains  only  white  noise.  The  purpose  of  the  prewhitening  filter  is 
to  remove  the  1/f  noise  generated  by  the  measurement  equipment.  To  this  end,  we  must 
first  obtain  an  estimate  of  the  noise  and  then  use  that  estimate  to  generate  an  appropriate 
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filter.  Since  this  detection  algorithm  sweeps  a  single  tone  across  a  bandwidth  of  interest,  it 
is  assumed  that  the  prewhitening  filter  is  only  acting  on  one  sinusoidal  frequency  embedded 
in  pink  noise  at  any  given  time.  It  has  been  shown  though,  that  this  filter  is  capable  of 
whitening  any  noise  spectrum  containing  narrowband  signals.  Proceeding  with  a  DFT  of 
the  zero  mean  input  as  described  by  Equation  (2.31),  yields  an  output  spectrum  whose 
exponential  decreasing  shape  is  of  interest.  Recall  that  we  also  have  to  preserve  the  single 
sinusoid  spectrum  (spike)  as  best  we  can.  To  this  end,  a  sixth  order  least-square  polynomial 
approximation  is  adopted  as  an  estimator  and  has  a  solution  of  the  form  [24] 

Nestif)  =  do  +  aif  +  ■  •  •  +  akfk,  (2.36) 


where  k  is  the  approximation  order,  and  used  in  the  following  derivations.  The  residual  of 
Equation  (2.36)  is  given  by 

N 

R2  =  y~^[iV(qc  —  (do  +  aifi  +  •  •  •  +  d&/f)]2,  (2.37) 

i=  1 

where  N^\c  is  the  ith  DFT  noise  coefficient.  To  minimize  the  residual,  its  partial  derivatives 
with  respect  to  each  of  the  approximation  coefficients  are  obtained  and  set  to  zero.  This 
yields  a  Vandermonde  matrix  relating  these  coefficients  to  the  actual  noise  values.  These 
partial  derivatives  are  expressed  as 

* 


d(R ' 


8a0 

d{R2 


-  —  —2  5>Wc  —  (ao  +  di/j  +  •  •  •  +  akfi)]  —  0 


i=  1 
N 


da\ 


~  —  ~ 2  y^[-^V(i)c  —  (do  +  di/j  +  . . .  +  d kft)]f  —  0 


i=  1 


d(R2 

dai. 


N 


■2  y^[-V(qc  —  (ao  +  di  fi  +  ...  +  akfi)]fk  —  0. 


(2.38) 


i= 1 


taking  the  derivatives  in  Equation  (2.38)  and  distributing  the  summations  leads  to 

N  N  N 

aon  +  di  ^  fi  +  . . .  +  dfc  ^  /jA  =  ^  N(i)c 


i= 1  i= 1 

N  N 


i= 1 

N  N 

ao  ^  fi  +  dl  ^  f2  +  ...  +  dfc  ^  /A+1  =  ^  fiN(i)c 
i=  1  1=1  i= 1  i=  1 

N  N  N  N 

a0J2ft  +  aiJ2fi+1  +  ---+^J2fik  =  J2fiN(i)c- 

i=  1  i=l  i= 1  i= 1 


(2.39) 
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Equation  (2.39)  can  be  written  in  matrix  form  as 

£2=1  fi  ■■■  £,= if?  ]  \  «o  1  [  EhN[c)i 

EN  r  r2  rk+ 1  at  n 

i=  1  Ji  2-ji=\  Ji  *  *  *  Z-a=l  Ji  2-^%—  1  (c)i^  ^2 

rk  sr^N  rk+1  sr^N  r2k  n ,  /y  fk 

L  2^i=l  Ji  2^=1  Ji  2^i=lJi  J  l  ak  J  L  Z^i=l  iy,(c)iJi  J 

Simplifying  the  notation,  we  can  obtain  the  matrix  for  a  least  squares  fit  by  writing 

i  /i  •  •  •  /f  ]  Uo  i  r  iv(C)i 

1  /2  •••  /2fe  ax  _  iV(c)2 

.  1  /n  •  •  •  fn  .  .  ak  _  .  W(c)n  J 

Using  matrix  notation  to  represent  Equation  (2.41),  the  equation  for  a  polynomial  fit  is 
given  by 

y  =  la.  (2.42) 

The  solution  results  by  premultiplying  both  sides  by  the  transpose  of  X 

XTy  =  XTX&,  (2.43) 

to  yield  a,  as 

a  =  ( XTX)~1XTy .  (2.44) 

A  sixth  order  least  squares  polynomial  approximation  was  heuristically  chosen  to 
provide  a  close  fit  to  the  noise  curve  while  minimizing  the  order  to  not  follow  the  sinusoidal 
spike  that  is  being  swept  from  0-10  kHz.  Once  an  approximation  of  the  colored  noise  is 
achieved,  the  reciprocal  is  taken  and  normalized  to  have  maximum  gain  of  one.  This  is 
mathematically  expressed  as 

fj  m  HM) 
wU)  Max{Hw(f)} 

Hw(f )  =  Nest(f)-\  (2.45) 

where  Hw(f)  is  the  normalized  fitler.  Finally,  the  filter  is  applied  to  the  input  signal 
producing  a  sinusoid  in  white  noise  of  slightly  decreased  amplitude  as  can  be  seen  in  Figure 
2.6. 


(2.41) 
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Frequency  (Hz)  Frequency  (Hz) 


Figure  2.7:  A  2  kHz  sinusoid  in  colored  pink  noise  (A),  prewhitening  filter  (B),  whitened 
sinusoid  (C),  and  whitened  sinusoid  with  corrected  amplitude  (D). 

2.2.6  Amplitude  Correction 

As  can  be  seen  in  Figure  2.6,  the  prewhitening  filter  attenuates  the  amplitude  of  the 
signal  being  filtered.  For  pink  noise,  the  algorithm  generates  a  high-pass  prewhitening  filter 
that  has  very  gentle  slope  in  what  amounts  to  the  transition  band.  If  the  signal  of  interest 
falls  in  this  area,  its  amplitude  will  be  greatly  decreased  and  will  not  be  accurately  reflected 
in  the  output  of  the  MUSIC  analysis.  To  compensate  for  this,  a  correction  coefficient  is 
introduced  to  rescale  the  amplitude  of  the  sinusoid  being  observed. 

Figure  2.5  depicts  how  the  correction  coefficient  is  generated.  The  magnitude  of 
the  sinusoid  in  the  original  DFT  is  divided  by  the  magnitude  of  the  sinusoid  in  the  whitened 
DFT.  This  ratio  is  then  multiplied  by  the  magnitude  of  the  sinusoid  in  the  whitened  spec¬ 
trum  to  return  it  to  the  original  value.  If  the  sinusoid  being  investigated  is  not  greater  than 
the  noise  floor,  the  prewhitening  assumes  its  magnitude  is  that  of  the  noise  floor,  and  still 
produces  a  sinusoid  in  white  noise  of  that  amplitude  and  at  that  frequency. 
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2.3  Multiple  Signal  Classification  (MUSIC)  Algorithm 

The  MUSIC  algorithm  was  initially  developed  by  Schmidt  [25],  Bienvenu,  and 
Kopp[26]  to  detect  emitter  location  by  direction  of  arrival  (DOA)  in  sensor  arrays.  Schmidt 
proved  that  the  MUSIC  algorithm  was  asymptotically  unbiased  in  estimating  1)  the  number 
of  incident  wavefronts;  2)  DOA;  3)  strength  and  cross  correlation  of  incident  waveforms; 
4)  noise/interference  ratio.  As  a  bonus,  the  MUSIC  algorithm  is  also  useful  as  a  frequency 
estimator  on  time  series  data[25].  Furthermore,  it  has  the  ability  to  provide  an  unbiased 
estimate  of  the  null-spectrum  of  uncorrelated  data  irregardless  of  the  SNR[27],  thus  making 
it  very  tolerant  to  noisy  applications. 

The  question  might  arise  as  to  why  the  MUSIC  algorithm  is  being  chosen  over  other 
spectral  methods  such  as  the  DFT,  maximum  likelihood  estimate  (MLE),  or  minimum  mean 
square  error  (MMSE).  The  DFT  is  the  most  well  known  for  its  ease  of  use,  but  also  for 
its  shortcomings.  In  this  application,  it  is  the  necessity  to  capture  enormous  data  sets  to 
minimize  spectral  leakage  and  to  lower  the  noise  floor  that  disqualifies  the  DFT.  Figure  (2.7) 
depicts  just  how  tolerant  the  MUSIC  algorithm  is  to  noise  as  it  is  capable  of  estimating 
the  frequency  when  the  signal  seem  to  be  totally  corrupted.  MMSE  has  been  discussed 

(A)  DFT  (B)  MUSIC  Algorithm 


Figure  2.8:  Frequency  spectrum  of  1  kHz  signal  in  noise  computed  with  Discrete  Fourier 
Transform  (A)  and  MUSIC  algorithm  (B). 

in  the  literature  for  frequency  estimation,  but  closed  form  solutions  are  difficult  and  lead 
to  computationally  complex  algorithms [2 8].  As  a  result,  we  opt  for  the  MUSIC  algorithm, 
which  was  stated  by  Schmidt [25]  to  asymptotically  approach  the  Cramer- Rao  accuracy 
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bound  (CRB)  and  proven  by  Stoica[29]  to  outperform  MLE. 

To  understand  why  the  MUSIC  algorithm  was  chosen,  it  is  helpful  to  begin  with 
the  data  model[29] 

y(t)  =  A(0)x(t)  +e(t) 

A (6)  =  [a(wi),  a(u;2),  •  •  • ,  a(w„)].  (2.46) 

In  Equation  (2.46),  y(t)  G  Cmxl  is  the  noisy  data  vector,  x(t)  G  Cnxl  is  the  vector  of 
signal  amplitudes,  e(t)  G  Cmxl  is  the  additive  noise,  and  A(0)  G  Cmxn  is  a  Vandermonde 
matrix  of  signal  frequencies  where  9n  =  .  In  order  to  use  the  model  presented 

in  Equation  (2.46),  some  basic  assumptions  must  be  made.  Table  (2.1)  lists  the  assumptions 
that  are  used  for  both  MUSIC  and  MLE  algorithms.  In  order  to  estimate  the  frequencies 


Table  2.1:  Assumptions  needed  for  MUSIC  and  MLE  frequency  estimator  algorithms. 


Assumption 

MUSIC 

MLE 

1 

m  >  n,  and  the  vectors  a(co)  corresponding  to 
(n  + 1)  different  values  of  u>  are  linearly  indepen¬ 
dent 

Yes 

Yes 

2 

E{e-{t)}  =  0,  E{e(t)e*(t)}  =  crl,  and 

E{e(t)eT  (t)}  =  0 

Yes 

Yes 

3 

The  matrix  P  =  E{x(t)x*  (f)}  is  nonsingular 
(positive  definite),  and  N  >  m 

Yes 

No 

4 

E{e(t)e* (s)}  =  E{e{t)e1  (s)}  =  0  for  t  /  s,  and 
e(t)  is  Gaussian  distributed. 

No 

Yes 

of  N  complex  sine  waves  from  single  experiment  data,  we  can  use  the  model 

N 

Uk  =  ^  7 PejWpk  +  ek  k  =  1, . . . ,  n.  (2.47) 

p=  i 

Equation  (2.47)  may  be  expressed  as  the  data  model  of  the  MUSIC  algorithm  by  using  the 
following  notation,  and  letting  m  be  some  integer  greater  than  n  such  that 

y(t)  =  [yt  ■  ■■ym]T 
a(wt)  =  [lejUt 

x(t)  =  [7ie^lt...7JVe^t]T 
e(t)  =  [et...em]T 


t  =  1, . . . ,  n. 


(2.48) 
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In  the  single  experiment  case,  e(t)  and  e(s)  are  correlated  for  t  7^  s,  contradicting  assump¬ 
tion  4  and  making  the  MLE  algorithm  not  applicable.  Assuming  that  MLE  could  be  used 
to  estimate  the  frequencies,  Stoica  shows  that  MUSIC  still  reaches  the  CRB  faster  than 
MLE  for  data  sets  of  length  n,  with  m  sensors.  Looking  at  the  ratios  of  MUSIC  and  MLE 
to  the  CRB, 


varML(u)/varCR{u) ) 

1  +  m  *  SNR 

(2.49) 

varMu(u)/varCR(tu) 

(AM)"1 

SNR  ’ 

(2.50) 

we  see  that  each  is  dependent  on  the  signal  to  noise  ratio  and  m.  For  MLE  to  reach  the 
CRB,  both  m  and  n  must  approach  infinity,  which  is  clearly  unattainable.  However,  for 
MUSIC  to  reach  the  CRB,  m  must  only  be  greater  than  n  for  relatively  large  n.  It  is  known 
that  MUSIC  is  inefficient  for  correlated  signals,  but  in  this  case  of  single  observations,  this 
drawback  does  not  present  a  problem. 

Therefore,  we  are  left  with  the  MUSIC  algorithm,  where  assumptions  1,2,  and  3 
are  satisfied  by  the  model  above  as  long  as  2 N  <  n  +  1.  This  final  constraint  simply  states 
that  the  number  of  complex  sinusoids  one  is  trying  to  analyse  can  not  exceed  the  number 
of  measurement  samples,  which  is  easily  achieved. 

Finally,  we  proceed  with  analyzing  the  eigen- decomposition  of  the  data  covariance 
matrix,  which  is  the  essence  of  the  MUSIC  algorithm.  Assuming  that  conditions  1,2,  and  3 
are  met,  the  covariance  matrix  of  observation  vector  y(t)  is  given  by 

R  =  £{y(t)y*(t)}  =  A(0)PA*(0)  +  a2 1,  (2.51) 

where  P  =  E{x(t)x*(t)}  is  the  autocorrelation  of  the  measured  data,  I  is  an  [m  x  m] 
identity  matrix,  and  <r2  is  the  noise  variance.  Taking  the  eigen  decomposition  of  R,  we 
obtain  m  eigenvalues,  D,  and  an  [m  x  m]  matrix,  V ,  of  eigenvectors  which  span  the  column- 
space  of  A.  Taking  the  eigenvectors  corresponding  to  the  first  n  largest  eigenvalues,  we 
can  define  a  matrix,  Es,  that  spans  the  signal  subspace[30].  The  remaining  N  =  m  —  n 
eigenvectors  span  the  noise  subspace  and  are  defined  by  the  [m  x  N]  matrix  En[25].  Due 
to  the  orthogonality  of  the  frequency  vectors  a(w)  in  A (6)  and  the  noise  subspace,  we  can 
compute  the  spectrum  as  the  inverse  of  the  inner  product  given  by 

JM  = 


a*(w)ENE^a(u;) 


(2.52) 
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2.4  Summary 

We  now  dispose  of  the  necessary  tools  to  analyze  the  time  domain  data  captured 
in  our  experimental  problem.  We  began  with  a  data  model  that  employs  a  perturbation 
analysis  to  model  a  lossy,  nonlinear  plane  wave  propagating  in  the  positive  x  direction.  This 
model  has  been  shown  to  outperform  the  classical  Bessel-Fubini  approximation  for  n  >  3 
order  and  for  distances  beyond  x  =  xjj.  This  model  has  been  shown  to  work  in  simulation, 
but  anechoic  chamber  data  should  be  analyzed  to  confirm  these  solutions.  However,  the 
ERL  currently  lacks  equipment  at  this  time  capable  of  generating  high  amplitude  ultrasound 
at  the  correct  frequency  to  generate  this  phenomenon.  Second,  we  must  preprocess  the  data 
to  remove  1/f  “pink”  noise  introduced  by  the  measurement  equipment,  and  average  the 
data  so  that  it  has  zero  mean.  Finally,  it  has  been  shown  that  the  MUSIC  algorithm  is 
optimal  in  the  sense  that  it  efficiently  approaches  the  CRB  and  is  the  only  spectral  analysis 
tool  that  can  effectively  process  this  low  SNR  data. 
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Chapter  3 

Experimental  Setup 

3.1  Anechoic  Chamber 

It  is  interesting  to  note,  that  using  the  aforementioned  nonlinear  detection  tech¬ 
nique  requires  measurements,  which  in  turn  require  that  a  suitable  environment  be  con¬ 
structed.  Due  to  the  inefficient  demodulation  from  two  ultrasonic  frequencies  to  one  base¬ 
band  tone,  and  the  fact  that  this  difference  frequency  lies  right  in  the  middle  of  the  noisy 
sonic  spectrum,  the  Electronics  Research  Laboratory  (ERL)  opted  to  construct  an  anechoic 
chamber.  Given  our  rather  wider  interest  in  acoustic,  electromagnetic,  and  acoustically 
modulated  electromagnetic  phenomenon,  building  a  dual  purpose  chamber  was  concluded 
as  the  best  option.  The  purpose  of  this  chamber  is  to  1)  attenuate  acoustic  reflections  to 
better  emulate  free  space;  2)  attenuate  electromagnetic  reflections  to  better  emulate  free 
space;  3)  attenuate  both  acoustic  and  electromagnetic  transmissions  for  safety  reasons.  This 
final  objective  is  necessary  due  to  the  high  input  power  levels  necessary  to  generate  mea¬ 
surable  nonlinear  behavior.  For  example,  it  takes  120  dB  sound  pressure  level  (SPL)  of 
ultrasonic  energy  to  produce  66  dB  of  audible  sound  at  1  m.  Typically,  acoustic  energy 
at  this  amplitude  causes  instant  permanent  hearing  damage.  While  most  scientists  believe 
that  high  amplitude  ultrasound  is  harmless  to  humans,  there  is  little  research  confirming 
its  effects  in  confined  spaces  or  with  directed  sources.  Furthermore,  acoustic  transmission 
attenuation  is  a  byproduct  of  the  anechoic  chamber’s  RF  shielding,  which  is  necessary  due 
to  the  well  known  dangers  of  high  energy  microwaves. 

Since  ERL  is  interested  in  measuring  life-sized  targets,  the  anechoic  chamber  was 
constructed  as  large  as  practically  possible.  Overall  external  chamber  dimensions  are  96 
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inches  in  width,  72  inches  in  height,  and  144  inches  in  length.  Due  to  the  thickness  of 
absorbent  materials  used  in  construction,  the  internal  usable  dimensions  are  76  inches  in 
width,  52  inches  in  height,  and  120  inches  in  length.  Chamber  construction  began  by 
building  a  raised  floor  to  provide  wire  runs  to  various  test  equipment  located  within  the 
chamber.  The  floor  measures  109  inches  in  width  by  156  inches  in  length  and  rests  on  top 
of  9”x4”  posts.  Starting  from  the  lower  most  layer,  the  floor  is  comprised  of  one  layer 
3/4”  plywood,  one  layer  cement  board,  one  layer  copper  mesh,  and  two  layers  of  6.0  mm 
thick  Acoustiblok.  The  walls  and  ceiling  of  the  anechoic  chamber  are  all  constructed  using 
the  same  process.  Support  for  the  walls  is  by  way  of  an  extruded  aluminum  space  frame 
manufactured  by  80/20  Corporation.  The  outermost  layer  of  each  wall  and  ceiling  panel  is 
comprised  of  copper  mesh,  manufactured  by  TWP  Inc.  [31],  forming  a  Faraday  cage  around 
the  entire  chamber.  An  exploded  view  of  the  anechoic  chamber  can  be  seen  in  Figure  3.1. 

Forming  the  foundation  of  each  wall  and  ceiling  panel  is  a  layer  of  3.0  mm  thick 
Acoustiblok.  This  high-density  rubberized  material  provides  almost  2/3  of  the  through- wall 
attenuation  above  1.0  kHz  [32].  The  Acoustiblok  sandwiches  the  copper  mesh  against  the 
supporting  frame  and  is  held  in  place  with  1.5  in.  nylon  bolts  spaced  4  in.  apart.  At  every 
seam  in  the  Acoustiblok,  acoustical  sealant  and  tape  is  used  to  further  improve  soundproof¬ 
ing.  Attached  to  the  inside  surface  of  the  Acoustiblok  are  2  ft.  by  4  ft.  panels  of  QuietBoard 
glued  to  Melamine  foam,  both  manufactured  by  American  Micro  Industries.  Each  panel  is 
attached  by  six  nylon  bolts,  and  forms  the  surface  to  which  the  RF  absorbing  foam  tiles 
are  glued.  Acoustiblok  sealant  is  again  used  at  all  panel  joints  to  ensure  soundproofing. 
The  innermost  layer  consists  of  RF  absorbing  geometric  tiles.  Two  types  of  this  tile,  each 
measuring  2  ft.  square,  were  used.  The  pyramidal  style  tile  was  used  in  the  most  sensitive 
areas,  such  as  along  the  back  wall  where  the  most  intense  energy  will  accumulate.  Eggshell 
tiles  were  used  to  fill  in  the  remaining  spaces.  All  of  the  tiles  were  arranged  in  a  manner  that 
reduced  the  possibility  of  generating  standing  waves  under  continuous  signal  generation. 

Acoustic  performance  of  the  ERL  anechoic  chamber  is  based  on  transmission  and 
reflection  loss  given  by  [33] 

Lt  =  201og10(^  (3.1) 

Lr  =  2Ologlo0|V  (3.2) 

where  Pt,  Pr,  and  Pi  represents  the  acoustic  pressure  of  the  transmitted,  reflected,  and 
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Figure  3.1:  Exploded  view  of  chamber  frame  and  floor.  Magnified  cross  section  shows 
wall  construction  detail.  The  Acoustiblok  and  copper  mesh  are  held  in  place  with  0.25  x 
1.5  in.  nylon  flat  head  bolts  and  nuts  attaching  these  layers  to  each  angle  bracket.  The 
Melamine /Quiet  Board  panels  are  held  in  place  with  0.25  x  4  in.  nylon  carriage  bolts. 
Finally,  RF  tiles  are  glued  to  the  Quiet  Board  using  contact  cement. 


incident  waves  respectively.  These  parameters  are  equivalent  to  those  obtained  from  a  S- 
parameter  matrix  of  a  two  port  device.  For  these  equations,  a  factor  of  20  is  used  since 
pressure  is  the  electrical  equivalent  of  voltage,  and  power  is  voltage  squared.  Figure  3.2 
depicts  the  acoustic  performance  of  the  ERL  anechoic  chamber.  RF  characterization  was 
also  performed,  but  is  not  included  in  this  document  since  it  does  not  pertain  to  acoustics. 


3.2  Test  Setup 


Two  excitation  signals,  of  frequency  f\  and  fi  are  produced  using  two  Marconi 
2024  signal  generators,  and  are  amplified/transmitted  using  Audio  Spotlight  parametric 
arrays  from  Holosonic  as  depicted  in  Figure  3.3.  For  these  experiments,  the  reflected  wave 
microphone  and  parametric  arrays  are  located  approximately  2.0  meters  from  the  target. 
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Figure  3.2:  Acoustic  transmission  and  reflection  loss  of  back  wall  of  anechoic  chamber. 
Above  50  kHz  where  higher  power  signals  can  be  generated,  the  periodicity  of  the  attenua¬ 
tion  indicates  that  a  resonance  is  being  generated. 

The  stationary  tone  f±  is  kept  at  55  kHz  while  /2  is  swept  from  55-65  kHz  in  100  Hz 
increments. 

Measurements  of  the  incident  and  reflected  sound  pressure  amplitudes  are  recorded 
using  two  377B01  condenser  microphones  from  Piezotronics,  and  digitized  by  a  PXI-5922 
High  Speed  Digitizer.  One  microphone  is  used  to  record  the  incident  wave,  while  the  other 
measures  the  reflected  wave.  A  sample  rate  of  50  kHz  is  used  and  50,000  time  samples  are 
captured  for  each  frequency  increment.  A  summary  of  the  test  equipment  specifications 
and  part  numbers  may  be  found  in  Appendix  H. 

Targets  chosen  for  this  experiment  include  glass,  metal,  wood,  and  fiberglass  con¬ 
tainers  of  both  regular  and  irregular  shapes  as  shown  in  Figure  3.4.  Targets  of  regular  shape 
were  perfect  10-inch  cubes.  Both  the  glass  and  wooden  cubes  have  sidewall  thickness  of 
approximately  1/4-inch,  while  the  metal  cube  is  constructed  of  20-gauge  steel.  An  irregular 
shaped  target,  emulating  a  piece  of  stone  allowed  us  to  test  a  more  realistic  target.  It  is 
constructed  of  molded  fiberglass  measuring  1/16- inch  thick,  with  the  solid  variant  filled 
with  concrete  in  lieu  of  sand  stone  at  2.3  g/crn3. 

To  single  out  the  one  phenomenon  we  are  trying  to  capture,  solid  and  hollow 
versions  of  the  same  target  were  placed  in  identical  positions  on  the  test  stand.  Usually,  a 
hollow  target  can  be  filled  with  some  material  (sand  in  these  experiments)  without  moving 
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Figure  3.3:  Top  view  of  anechoic  chamber  showing  two  Audio  Spotlight  transmitters  and 
two  condenser  microphones.  The  incident  microphone  is  only  used  to  monitor  incident 
sound  pressure  levels  and  is  not  necessary  for  detection. 


it  on  the  test  stand.  This  ensures  that  it  is  the  target’s  composition  that  produces  a  given 
acoustic  signature,  and  not  a  slight  change  in  its  alignment. 

At  the  surface  of  the  test  object,  the  two  tones  mix  to  produce  inter-modulation 
products.  The  transmitted  signals  are  shown  in  Figure  3.5  and  may  be  modeled  using[34] 

A1(t)  =  |Ai|  cos(a>if  +  <j) i)  +  \X2\  cos (uj2t  +  <fi2),  (3.3) 

to  yield  second  order  inter-modulation  products  of 

X\t)=(^j  [X2ej2w1t  +  2X1X*1 
+2X1X2ej{wi+W2)t  +  2XlX^ej{wi~W2)t 
+(Xl)2e~j2wit  +  2XlX2e^W2~wl)t 
+2XlX*2e-j(-wl+W2)t  +  X$ej2w2t 

+2X2X2  +  {X^)2e~j2w2t],  (3.4) 

where  Xn  is  the  amplitude  of  the  original  cosine  at  /j  and  f2,  producing  new  frequencies 
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Figure  3.4:  Targets  used  include  10-inch  cubes  of  glass  (A),  metal  (B),  wood  (C),  and 
simulated  sandstone  (D). 

at  /2  —  fi  and  /2  +  /i  as  well  as  third  order  inter-modulation  tones  and  higher  frequency 
tones  which  are  not  of  interest.  The  base  band  tone  at  frequency  /2  —  fi  is  swept  from  0.1- 
10  kHz  to  provide  a  frequency  response  or  acoustic  signature  of  the  target.  This  difference 
frequency  may  also  couple  into  the  target  object  where  it  can  potentially  excite  a  resonance. 


3.3  Measurements  and  Signal  Processing 

Upon  appending  the  time  domain  measurements  at  each  frequency  increment  to 
the  same  LabView  data  file,  we  proceed  with  the  signal  processing  techniques  discussed  in 
Chapter  2  to  analyze  the  data.  Given  the  large  data  sets,  and  the  fact  that  at  this  stage  of 
research,  we  are  monitoring  four  channels  of  data,  a  Matlab  graphical  user  interface  (GUI) 
environment  was  developed  to  aid  in  the  signal  processing.  This  program  allows  the  user 
to  easily  monitor  all  four  signals,  and  to  change  the  parameters  used  to  implement  the 
prewhitening  filter  and  MUSIC  algorithm.  It  also  provides  real  time  graphical  algorithm 
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Figure  3.5:  Spectrum  showing  excitation  signals  /i  =  55  kHz,  fi  =  65  kHz,  and  difference 
frequency  f\  =  10  kHz  generated  in  air. 

feedback  so  the  user  may  evaluate  the  noise  model  and  terminate  the  program  as  needed, 
e.g.,  a  parameter  needs  to  be  changed.  Figure  3.6  shows  a  screenshot  of  the  “Music  Gui.rn” 
program  and  Tables  3.1  and  3.2  define  the  inputs  and  outputs  respectively.  The  code  used 
to  implement  “Music_Gui.m”  is  available  in  Appendices  D-G. 

Upon  entering  all  of  the  correct  parameters  and  input  files,  a  user  may  begin 
to  analyze  the  data  by  simply  initiating  the  “start”  button.  The  four  clusters  of  two 
graphs  show  the  noise  model  used  to  implement  the  prewhitening  filter  and  the  MUSIC 
coefficients,  J(w),  for  hollow-incident,  hollow-reflected,  solid-incident,  and  solid-reflected 
signals  respectively.  The  prewhitening  filter  graphs  show  the  noisy  data  in  red,  the  noise 
approximation  in  green,  and  the  filtered  data  in  blue.  The  MUSIC  coefficients  graph  shows 
the  total  MUSIC  pseudospectrum,  J(w),  in  blue  and  the  coefficients,  J(uin),  chosen  by 
the  algorithm  in  red.  Ideally,  the  coefficients  shown  as  a  red  stem  plot  should  match  the 
peak  in  the  MUSIC  pseudospectrum.  However,  if  the  algorithm  can  not  find  a  peak  in  the 
pseudospectrum  at  the  frequency  being  examined,  the  coefficient  is  replaced  by  the  noise 


42 


Figure  3.6:  “Music_GUI.m”  signal  processing  application. 

floor  of  the  pseudospectrum,  which  is  approximately  30.  When  this  occurs,  the  zeroed 
coefficients  counter  is  incremented  to  also  update  the  percentage  of  valid  coefficients.  This 
may  be  seen  in  the  last  MUSIC  graph  of  Figure  3.6,  where  the  MUSIC  coefficient  does  not 
match  the  pseudospectrum.  Finally,  the  individual  J(con)  are  then  combined  to  generate  a 
real  time  frequency  response  plot,  S(f),  located  below  the  input  fields. 

The  program  “MUSIC_GUI.m”  described  above  produces  an  acoustic  signature  of 
both  the  hollow  and  solid  variants  of  the  same  target.  By  inspection,  it  may  be  seen  whether 
these  targets  respond  differently  if  they  are  hollow  or  solid,  thus  showing  that  a  detection 
algorithm  could  be  implemented.  While  the  scope  of  this  thesis  is  not  to  derive  a  detection 
algorithm,  the  “MUSIC_GUI.m”  program  does  produce  an  output  file  so  that  further  signal 
processing  may  be  pursued.  This  may  include  plotting  the  data  using  different  scales,  or 
taking  the  integral  of  the  acoustic  signature,  S(f),  as  shown  in  Chapter  4,  to  determine  the 
average  value,  which  may  be  a  good  indicator  of  the  targets  composition. 
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Table  3.1:  MUSIC  GUI  input  parameter  definitions. 


Parameter 

Definition 

FI 

Stationary  frequency  f\. 

F2  Start 

f-2  frequency  sweep  start. 

F2  Increment 

f'2  frequency  sweep  increment. 

F2  Stop 

/2  frequency  sweep  stop. 

Sample  Rate 

Sample  rate  used  by  LabView  Virtual  Instrument 
(VI). 

Samples 

Samples  captured  per  frequency  increment  by  Lab- 
View  VI,  used  to  parse  input  data  file. 

Hollow  Target 

Filename  of  hollow  target  data  file  (must  specify  path 
name) . 

Solid  Target 

Filename  of  solid  target  data  file  (must  specify  path 
name) . 

Output 

Output  file  name  (saved  to  same  location  as  MU- 
SIC.GULm). 

LSE  Order 

Least  square  estimate  order  for  noise  model. 

Filter  Samples 

The  number  of  samples  used  by  DFT  to  approximate 
noise  (must  be  less  than  Samples). 

Filter  Start  Frequency 

Start  frequency  of  prewhitening  filter  (should  be  zero 
for  pink  noise). 

Filter  Stop  Frequency 

Stop  frequency  of  prewhitening  filter. 

Y-Max,  Y-Min 

Plot  parameters  for  viewing  prewhitening  filter  and 
MUSIC  output  (X-Min,  X-Max  defined  by  range  of 
frequency  sweep. 

Number  of  Sinusoids 

The  number  of  sinusoids  the  MUSIC  algorithm  is  look¬ 
ing  for. 

Segment  Length 

The  length  of  partitioned  blocks. 

Overlap  Percent 

The  percentage  of  overlap  between  consecutive  blocks. 

Threshold 

Specifies  the  cutoff  for  the  signal  and  noise  subspace 
separation. 

Table  3.2:  MUSIC  GUI  output  parameter  definitions. 


Parameter 

Definition 

Index 

The  current  loop  index. 

Frequency 

The  current  frequency  increment. 

CC 

The  current  correction  coefficient. 

z_c 

The  total  number  of  zeroed  MUSIC  coefficients. 

% 

The  percentage  of  valid  MUSIC  coefficients. 
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Chapter  4 

Results  and  Conclusions 

4.1  Results 

Three  out  of  the  five  targets  tested  in  the  course  of  this  research  effort  show 
a  significant  difference  between  hollow  and  solid  variants.  These  targets  include  a  metal 
cube,  glass  cube,  and  simulated  rock.  The  plastic  container  shown  in  Figure  4.5  shows  some 
difference  in  solid  versus  hollow  variants,  while  the  wood  cube  shown  in  Figure  4.3  shows 
practically  no  difference.  To  verify  the  performance  of  the  MUSIC  algorithm,  a  constant 
amplitude  frequency  sweep  is  applied.  The  expected  flat  spectrum,  shown  in  Figure  4.1, 
validates  the  algorithm’s  performance.  Most  of  these  acoustic  signatures  were  generated 
with  only  2,000  MUSIC  samples,  as  opposed  to  the  50,000  samples  that  were  captured  and 
used  to  generate  similar  results  using  only  the  DFT.  This  shows  that  MUSIC  analysis  can 
produce  equivalent  acoustic  spectra  using  significantly  less  data  and  in  far  less  time.  Table 
4.1  summarizes  these  results  and  provides  a  metric  for  the  difference  between  hollow  and 
solid  variants.  This  metric  is  given  as  a  percentage  and  is  computed  by 

difference  =  hr  x  100%,  (4.1) 

where  Scut0ff  is  the  noise  floor  of  the  MUSIC  pseudospectrum.  In  Equation  4.1,  Sfir ,  Ssr , 
and  Si  are  the  frequency  response  average  values  of  the  hollow  reflected,  solid  reflected,  and 
incident  frequency  responses  respectively.  The  average  value  of  each  response  is  computed 
by 

1  rfstop 

Sav9  =  7 - —7 -  /  S(f)df,  (4.2) 

J  stop  J start  J f start 
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where  S(f)  is  the  frequency  response  of  the  target. 

Using  the  above  metric,  the  metal  cube  displays  the  greatest  difference  between 
solid  and  hollow  variants,  as  may  be  seen  in  Figure  4.2.  Figure  4.2  shows  that  the  empty 
metal  cube  reflects  less  acoustic  energy.  This  is  in  contrast  to  what  appears  in  the  responses 
of  the  glass  cube,  simulated  rock,  and  plastic  containers,  whose  empty  responses  show 
that  more  acoustic  energy  is  reflected.  Material  properties  that  can  affect  the  amount  of 
reflected  acoustic  energy  include  elasticity,  acoustic  impedance,  and  homogeneity  of  the 
target.  The  difference  between  a  metal  cube  and  other  targets  is  most  likely  based  on 
their  structural  make.  The  glass  cube,  simulated  rock,  and  plastic  container  are  all  molded 
from  one  continuous,  homogenous  material.  This  means  that  acoustic  energy  incident  on 
these  hollow  targets  gets  trapped  in  the  shell  and  can  not  pass  into  the  hollow  interior  due 
to  the  large  impedance  mismatch  between  solids  and  air.  On  solid  targets  of  this  type, 
the  mechanical  interaction  between  the  shell  and  interior  medium  allows  acoustic  energy  to 
propagate  into  the  interior  where  it  is  dissipated.  The  wood  and  metal  cubes  are  constructed 
from  individual  pieces  held  together  with  mechanical  fasteners.  This  may  prevent  trapped 
absorbed  acoustic  energy  in  the  shell  of  a  hollow  target  from  traveling  as  easily  around 
the  perimeter.  Furthermore,  the  metal  cube  has  significantly  greater  acoustic  impedance 
as  compared  to  the  other  targets  while  still  being  compliant.  This  is  in  contrast  to  the 
glass  cube  and  simulated  rock  which  have  low  impedance  and  high  rigidity.  The  material 
properties  are  also  summarized  in  Table  4.1. 
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Figure  4.1:  MUSIC  analysis  of  con¬ 
stant  amplitude  frequency  sweep  to 


test  algorithm. 


Metal  Cube 


Figure  4.2:  Metal  cube  shows  the 
greatest  difference  in  hollow  and 
solid  variants. 
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Wood  Cube 


Figure  4.3:  Wood  cube  shows  no  dif¬ 
ference  in  hollow  and  solid  variants. 


Glass  Cube 


Figure  4.4:  Glass  cube  shows  differ¬ 
ence  in  solid  and  hollow  variants,  but 
required  more  samples  to  do  so. 


One  important  observation  to  make  about  Figure  4.6  is  the  existence  of  a  resonant 
behavior  in  a  hollow  target.  These  are  evident  in  the  equally  spaced  dips  in  the  reflected 
sound  energy.  In  Figure  4.6,  these  dips  fall  at  intervals  of  3.5  kHz.  Knowing  the  speed  of 
sound  in  the  fiberglass  shell  to  be  2,650  m/s,  one  can  calculate  the  resonant  length  to  be 


t  = 

Vr 


(4.3) 


Here  vs  is  the  sound  speed  in  the  material,  and  fr  is  the  resonant  frequency.  For  the  hollow 
rock,  Equation  (4.3)  yields  a  resonant  length  of  0.37  nr,  which  is  very  close  to  the  actual 
rock  size  of  0.40  nr.  If  such  a  resonant  response  can  be  found  in  all  targets,  it  will  be  possible 
to  estimate  not  only  the  composition  of  a  target,  but  also  its  size. 


4.2  Future  Work 

While  the  results  for  this  research  are  thus  far  promising,  there  is  still  much  work 
to  be  done  in  the  way  of  improved  sensing,  increased  SNR,  and  expanded  functionality. 
First  of  all,  the  distinction  should  be  made  between  mechanical  and  acoustic  resonance. 
Mechanical  resonance  is  the  movement  of  the  target’s  surface  due  to  some  applied  harmonic 
force  and  is  given  by  Equation  (4.3).  Mechanical  resonance  is  difficult  to  measure  with 
microphones  due  to  the  poor  coupling  between  solids  and  air,  and  could  be  more  easily 
measured  with  equipment  such  as  a  Laser  Doppler  Vibrometer  (LDV).  Acoustic  resonance 
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Plastic  Container 


Figure  4.5:  Plastic  container  shows 
some  difference  in  hollow  and  solid 
variants. 


Simulated  Rock 


Figure  4.6:  Simulated  rock  shows 
significant  difference  in  hollow  and 
solid  variants  and  possible  resonant 
response. 


is  the  movement  of  air  molecules  within  a  target,  which  if  modeled  as  a  Helmholtz  resonator 
has  a  resonant  frequency  given  by  [35] 

fr=2n{vi  (44) 

In  Equation  (4.4),  c,  S,  l,  and  V  are  the  sound  velocity,  cross  sectional  area  of  the  neck, 
length  of  the  neck,  and  volume  of  the  cavity,  respectively.  Since  the  resonant  frequency  is 
inversely  proportional  to  the  volume,  acoustic  resonance  can  give  more  precise  indication 
of  the  amount  of  substance  in  a  target  in  addition  to  it  being  hollow  or  solid. 

Another  improvement  that  can  be  made  is  to  increase  the  transmit  power  of  the 
ultrasonic  sources.  The  amplifiers  and  transducers  from  Holosonic  used  in  this  research  have 
a  limit  of  13.6  Vpp  input.  This  limit  is  primarily  due  to  the  large  operating  bandwidth  of  10 
kHz.  If  smaller  bandwidth  transducers  are  used,  such  as  those  from  Airmar,  input  signals 
as  high  as  1000  Vpp  can  be  used.  These  transducers  have  a  much  narrower  bandwidth 
of  approximately  four  percent  of  the  center  frequency,  and  must  be  operated  at  low  duty 
cycles.  Furthermore,  the  use  of  pulsed  waveforms  such  as  radar  can  increase  the  relative  in¬ 
stantaneous  power  applied  to  the  target.  Linear  frequency  modulated  (LFM)  chirps  having 
sufficiently  wide  bandwidth  in  the  frequency  domain,  simulate  a  delta  function  in  the  time 
domain  providing  a  more  realistic  impulse  response  of  a  target. 

In  addition,  using  an  LFM  chirp  as  in  acoustic  radar  can  provide  range  data 
as  well  as  target  composition.  Using  the  well  known  method  of  stretch  processing  (pulse 
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Table  4.1:  Summary  of  results  for  hollow  targets  (HT)  and  solid  targets  (ST). 


Target 

MUSIC 

Samples 

Elasticity 

(GPa) 

Acoustic 

Impedance 

(MRayl) 

Avg. 

HT 

Amp. 

Avg. 

ST 

Amp. 

%  Diff. 

Simulation 

2000 

N/A 

N/A 

36.61 

36.61 

N/A 

Metal  Cube 

2000 

190-210 

47 

30.46 

33.89 

53.37 

Simulated 

Rock 

2000 

45 

2.85 

34.98 

31.86 

49.68 

Glass  Cube 

9000 

65-90 

12.5 

34.40 

32.02 

35.87 

Plastic 

Container 

2000 

0.8 

2.33 

32.78 

31.90 

13.87 

Wood  Cube 

9000 

8-13 

1. 5-3.0 

34.53 

34.54 

0.24 

compression),  a  high  resolution  range  spectrum  can  be  produced  by  comparing  the  reflected 
LFM  chirp  to  the  transmitted  one.  In  stretch  processing,  the  transmitted  signal  is  modeled 


as  [36] 


St(t)  =  COS  (27 r  (f0t  +  7^2))  , 


(4.5) 


where  fo  is  the  LFM  start  frequency,  and  n  =  Bjr  is  the  LFM  bandwidth  divided  by  the 
pulse  duration.  If  we  introduce  a  time  delay  At  =  2 R/c  and  attenuation  coefficient  a  for  a 
radar  cross  section  (RCS),  antenna  gain,  and  range  attenuation,  the  received  signal  is  given 
by 


sr(t)  =  a  cos 


2vr  (f0{t  -  At)  +  |(t  -  At)2) 


(4.6) 


Mixing  of  the  transmitted  and  received  waveforms  produces  a  signal  with  an  instantaneous 
frequency  that  is  dependent  on  the  range(s)  of  the  target.  Assuming  a  peak  in  the  frequency 
response  at  f  \ ,  the  range  of  the  target  is  computed  by 


Ri  = 


(4.7) 


Using  the  specifications  of  the  LabView  equipment  already  available  in  the  lab,  a  simulation 
of  a  stretch  processor  implemented  in  Matlab  shows  an  achievable  range  resolution  of  5  cm 
as  shown  in  Figure  4.7. 

In  summary,  the  use  of  chirped  signals  can  increase  the  effective  power  incident 
on  the  target,  while  providing  range  data.  With  the  addition  of  a  LDV,  mechanical  reso¬ 
nance  can  be  more  easily  measured,  leaving  the  standoff  microphones  to  focus  on  measuring 
acoustic  resonance,  which  provides  greater  volume  resolution. 
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Range  (m) 

Figure  4.7:  DFT  of  pulsed  compressed  received  signal  showing  resolvable  targets  at  1.5  and 
1.55  m. 

4.3  Conclusion 

A  novel  approach  using  nonlinear  ultrasonic  parametric  arrays  for  standoff  analysis 
of  targets  has  been  presented.  It  has  been  shown  that  this  approach  is  capable  of  distin¬ 
guishing  between  solid  and  hollow  variants  of  the  same  target,  with  the  potential  to  gain 
even  more  information  about  the  target  with  minimal  additional  complexity.  Furthermore, 
the  third  order  nonlinear  air  model  presented  in  Chapter  2,  performs  better  than  conven¬ 
tional  models  in  predicting  third  order  harmonic  growth  in  simulation.  This  model  will 
help  to  better  predict  nonlinear  ultrasonic  propagation  so  that  estimates  of  incident  energy 
on  targets  at  standoff  ranges  will  be  more  accurate.  The  combination  of  this  work,  with 
that  of  finite  element  models  investigated  by  Vetreno,  may  allow  faster,  higher  dimensional 
models  to  be  generated[37]. 

Furthermore,  using  the  MUSIC  algorithm  for  spectral  estimation,  provides  acoustic 
signatures  using  only  2,000  time  samples  as  opposed  to  50,000  for  Fourier  analysis.  This 
results  in  a  significant  decrease  in  computation  time  while  being  more  tolerant  to  low  SNR. 
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As  a  final  contribution,  the  “Air_ModeLGUI.m”  and  “MUSIC_GUI.m”  applications  have 
been  included  to  aid  in  future  development  and  provide  easier  algorithm  implementation. 
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Appendix  A.  Discrete  Fourier  Transform  Matlab  Code 

//  function  [Xk_mag,  frequencies,  k,  Xk_real,  Xk_imag]  =  general_ft (data, 
//  delta_t,  start_freq,  end_freq) 

// 


// "/ 

'0  Discrete 

Fourier  Transform 

// "/ 

f 

0 

// "/ 

usage : 

general_f t (S_t ,  delta_t,  start_freq,  end_freq) 

// "/ 

f 

0 

// "/ 

'0  inputs: 

S_t  =  time  domain  data 

// "/ 

f 

0 

delta_t  =  time  between  data  samples 

// "/ 

f 

0 

start_freq  =  start  of  Fourier  spectrum 

// "/ 

f 

0 

end_freq  =  end  of  Fourier  spectrum 

// "/ 

t 

0 

// "/ 

'0  outputs: 

:  Xk_mag  =  magnitude  spectrum 

// "/ 

t 

0 

frequencies  =  frequency  vector 

// "/ 

f 

0 

k  =  k  vector  (0  1  2  ...  (end_freq-start_freq)/df ) 

// "/ 

f 

0 

Xk_real  =  real  magnitude  spectrum 

// "/ 

/ 

0 

Xk_imag  =  imaginary  magnitude  spectrum 

// 

//  "/length  of  input  time  domain  data  signal 
//  N  =  length (data) ; 

//  "/(frequency  resolution  from  time/bandwidth  product 
//  df  =  1/ (N*delta_t) ; 

//  %total  time  duration  of  input  signal 
//  total_time  =  N*delta_t; 

// 

//  for  i=l : ( ( (end_f req-start_freq) /df )+l) 

//  "/generation  of  output  frequency  vector 

//  frequencies (i)  =  start_freq+((i-l)*df) ; 

//  "/generation  of  k  vector 

//  k(i)  =  frequencies (i)  *  total_time; 

//  Xk_real(i)  =  0; 

//  Xk_imag(i)  =  0; 

//  for  j=l:N 

//  Xk_real(i)  =  (Xk_real(i)+(data(j)*cos(-2 

//  *3 . 14159265358979323846* (j-l)*k(i)/N))) ; 

//  Xk_imag(i)  =  (Xk_imag(i)+(data(j)*sin(-2 

//  *3 . 14159265358979323846* (j-l)*k(i)/N))) ; 

//  end 

//  "/real  amplitude  correction 

//  Xk_real(i)  =  2*Xk_real(i)/N; 

//  "/imaginary  amplitude  correction 

//  Xk_imag(i)  =  2*Xk_imag(i)/N; 


//  "/.magnitude  spectrum 

//  Xk_mag(i)  =  sqrt  ( (Xk_real(i)  ~2)+(Xk_imag(i)  ~2) ) ; 

//  "/.conversion  to  dB  scale 

//  Xk_dB(i)  =  20*logl0((Xk_mag(i)*(sqrt(2)/2))/0. 00002) 

//  end 

//"/.uncomment  to  output  amplitude  in  dB  scale. 

//  0/„Xk_mag  =  Xk_dB ; 

#include  <math.h> 

#include  <stdio.h> 

#include  "mex.h" 

#include  "matrix. h" 

#include  <stdlib.h> 

void  mexFunction(int  nlhs,  mxArray  *plhs  []  , 
int  nrhs ,  const  mxArray  *prhs [] ) { 

double  *data; 
double  *d_t ; 
double  *f_st; 
double  *f_sp; 
double  *real; 
double  *imag; 
double  *mag; 
double  *freq; 
double  *k; 

int  i ,  j  ; 
int  N; 

double  total_time; 
double  samples; 
double  d_f ; 

double  Pi  =  3.14159265358979323846; 

data  =  mxGetPr (prhs  [0] ) ; 
d_t  =  mxGetPr (prhs [1] ) ; 
f_st  =  mxGetPr (prhs  [2] ) ; 
f_sp  =  mxGetPr (prhs [3] ) ; 

N  =  mxGetN(prhs [0] ) ; 
total_time  =  N*(*d_t); 


if  (*f _sp>*f _st){ 

d_f  =  l/total_time ; 
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samples  =  ((*f_sp-(*f_st))/d_f)+l; 

} 

else{ 

d_f  =  0; 
samples  =  1; 

> 

plhs [0]  =  mxCreateDoubleMatrix(l , samples ,mxREAL) ; 
plhs[l]  =  mxCreateDoubleMatrix(l , samples ,mxREAL) ; 
plhs  [2]  =  mxCreateDoubleMatrix(l , samples ,mxREAL) ; 
plhs  [3]  =  mxCreateDoubleMatrix(l , samples ,mxREAL) ; 
plhs  [4]  =  mxCreateDoubleMatrix(l , samples ,mxREAL) ; 

mag  =  mxGetPr (plhs [0] ) ; 
freq  =  mxGetPr (plhs [1] ) ; 
k  =  mxGetPr (plhs [2] ) ; 
real  =  mxGetPr (plhs [3] ) ; 
imag  =  mxGetPr (plhs [4] ) ; 

for(i=0;  i<samples;  i++){ 

freq[i]  =  (*f _st)+(i*d_f ) ; 
k[i]  =  freq[i] *total_time; 
real[i]  =  0; 
imag  [i]  =  0; 
f or (j=0;  j<N ;  j++){ 

real [i]  =  (real [i] +  (data[j] *cos (-2*Pi*j*k  [i] /N) ) ) ; 
imag[i]  =  (imag [i] +(data[j] *sin(-2*Pi*j *k [i] /N) ) ) ; 

> 

real  [i]  =  2*real[i]/N; 
imag[i]  =  2*imag[i]/N; 

mag[i]  =  sqrt(pow(real [i] ,2)+pow(imag[i]  ,2)) ; 
//"/oUncomment  to  output  magnitude  in  dB  scale 
//mag[i]  =  20*logl0((mag[i] *(sqrt(2)/2))/0. 00002) ; 

> 

> 
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Appendix  B.  Inverse  Discrete  Fourier  Transform  Matlab  Code 

//  function  [S_rec]  =  general_ift(Xk_real,Xk_imag,k,N) 

// 

//  %  Inverse  Discrete  Fourier  Transform 

//  7„ 

//  %  usage:  general_if t (Xk_real ,  Xk_imag,  k,  N) 

//  7» 

//  7o  inputs:  Xk_real  =  real  fourier  coefficients 

//  7o  Xk_imag  =  imaginary  fourier  coefficients 

//  7o  k  =  k  vector  (0  1  2  ...  (end_freq-start_freq)/df ) 

//  7o  N  =  number  of  output  time  samples 

//  7» 

//  7o  outputs:  S_rec  =  reconstructed  signal  of  length  N 

// 

//  for  n=l:N 
//  S_rec(n)  =  0; 

//  for  k_ind  =  l:length(k) 

//  S_rec(n)  =  S_rec(n)+((Xk_real(k_ind)*cos(2*pi*k(k_ind)* 

//  (n-1) /N) ) -(Xk_imag(k_ind) *sin(2*pi*k(k_ind) * (n-1) /N) ) ) ; 

//  end 

//  end 

#include  <math.h> 

#include  <stdio.h> 

#include  "mex.h" 

#include  "matrix. h" 

#include  <stdlib.h> 

void  mexFunction(int  nlhs,  mxArray  *plhs [] , 
int  nrhs ,  const  mxArray  *prhs [] ) { 

double  *Xk_real; 
double  *Xk_imag; 
double  *k; 
double  *N ; 
double  *Srec; 

int  i ,  j  ; 
int  samps; 

double  Pi  =  3.14159265358979323846; 

Xk_real  =  mxGetPr (prhs  [0] ) ; 

Xk_imag  =  mxGetPr (prhs  [1] ) ; 
k  =  mxGetPr (prhs [2] ) ; 
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N  =  mxGetPr (prhs [3] ) ; 
samps  =  mxGetN(prhs [0] ) ; 

plhs  [0]  =  mxCreateDoubleMatrix(l , *N,mxREAL) ; 

Srec  =  mxGetPr (plhs [0] ) ; 

for(i=0;  i<*N;  i++){ 

Srec  [i]  =  0; 

for(j=0;  j<samps;  j++){ 

Srec[i]  =  Srec [i]  +  ( (Xk_real [j] *cos (2*Pi*k  [j] *i/ (*N) ) ) 
- (Xk_imag [j]*sin(2*Pi*k[j]*i/(*N)))) ; 

> 

} 

} 
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Appendix  C.  Air  Model  GUI  Application  Matlab  Code 


function  varargout  =  Air_Model_GUI (varargin) 

7.  AIR_MODEL_GUI  M-file  for  Air_Model_GUI . f ig 

7.  AIR_MODEL_GUI ,  by  itself,  creates  a  new  AIR_MODEL_GUI  or  raises 

7o  the  existing  singleton*. 

1 

7o  H  =  AIR_MODEL_GUI  returns  the  handle  to  a  new  AIR_MODEL_GUI  or 

7.  the  handle  to  the  existing  singleton*. 

7. 

*/.  AIR_MODEL_GUI(’ CALLBACK’  ,hObject,  eventData,  handles,  ... )  calls  the 

%  local  function  named  CALLBACK  in  AIR_MODEL_GUI .M  with  the  given 

7o  input  arguments. 

7. 

%  AIR_MODEL_GUI ( ’Property’ , ’Value’ , . . .)  creates  a  new  AIR_MODEL_GUI 

7o  or  raises  the  existing  singleton*.  Starting  from  the  left, 

7o  property  value  pairs  are  applied  to  the  GUI  before 

7o  Air_Model_GUI_OpeningFunction  gets  called.  An  unrecognized 

7o  property  name  or  invalid  value  makes  property  application  stop. 

7o  All  inputs  are  passed  to  Air_Model_GUI_OpeningFcn  via  varargin. 
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l 

7.  *See  GUI  Options  on  GUIDE'S  Tools  menu.  Choose  "GUI  allows  only 

*/.  one  instance  to  run  (singleton) " . 

7. 

7.  See  also:  GUIDE,  GUIDATA,  GUIHANDLES 

7o  Edit  the  above  text  to  modify  the  response  to  help  Air_Model_GUI 

7.  Last  Modified  by  GUIDE  v2.5  25-Feb-2008  14:54:33 

7o  Begin  initialization  code  -  DO  NOT  EDIT 
gui_Singleton  =  1; 

gui_State  =  struct( ,gui_Name> ,  mfilename,  ... 

’ gui_Singleton’ ,  gui_Singleton,  . . . 

’ gui_OpeningFcn’ ,  @Air_Model_GUI_OpeningFcn,  ... 

’ gui_OutputFcn’ ,  @Air_Model_GUI_OutputFcn,  . . . 

’ gui_LayoutFcn’ ,  [],... 

’ gui_Callback'  ,  [] )  ; 

if  nargin  &&  ischar (varargin{l]-) 

gui_State  ,gui_Callback  =  str2func  (varargin{l]-) ; 

end 

if  nargout 

[varargout{l :nargout}]  =  gui_mainf cn(gui_State ,  varargin{ : }) ; 

else 

gui_mainf cn(gui_State ,  varargin{ : }) ; 

end 

7o  End  initialization  code  -  DO  NOT  EDIT 


7o - Executes  just  before  Air_Model_GUI  is  made  visible. 

function  Air_Model_GUI_OpeningFcn(hObject ,  eventdata,  handles,  varargin) 
7o  This  function  has  no  output  args,  see  OutputFcn. 

7o  hObject  handle  to  figure 

7o  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
7o  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

7o  varargin  command  line  arguments  to  Air_Model_GUI  (see  VARARGIN) 


7o  Choose  default  command  line  output  for  Air_Model_GUI 
handles . output  =  hObject; 


7o  Update  handles  structure 
guidata (hObject ,  handles); 
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%  UIWAIT  makes  Air_Model_GUI  wait  for  user  response  (see  UIRESUME) 
%  uiwait (handles . figurel) ; 


l  -  Outputs  from  this  function  are  returned  to  the  command  line. 

function  varargout  =  Air_Model_GUI_OutputFcn(hObject ,  eventdata,  handles) 
'/„  varargout  cell  array  for  returning  output  args  (see  VARARGOUT) ; 

*/.  hObject  handle  to  figure 

*/.  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

*/.  Get  default  command  line  output  from  handles  structure 
varargout {1}  =  handles . output ; 

°/„Create  SPL  Input 

function  editl_Callback(hObject ,  eventdata,  handles) 
handles. SPL  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

function  editl_CreateFcn(hObject ,  eventdata,  handles) 
handles. SPL  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , 'BackgroundColor ') , 
get (0, ' def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ' , ’ white ’ ) ; 

end 

°/„Create  frequency  input 

function  edit2_Callback (hObject ,  eventdata,  handles) 
handles. f  =  str2double (get (hObject ,  ’ String' )) ; 
guidata(hObject .handles) ; 

function  edit2_CreateFcn(h0bject ,  eventdata,  handles) 
handles. f  =  str2double (get (hObject ,  ’ String' )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , 'BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  Gamma  input 

function  edit3_Callback (hObject ,  eventdata,  handles) 
handles. Y  =  str2double (get (hObject String' )) ; 
guidata(hObject .handles) ; 

function  edit3_CreateFcn(h0bject ,  eventdata,  handles) 
handles. Y  =  str2double (get (hObject String' )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , 'BackgroundColor ’) , 
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get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

°/„Create  ideal  gas  constant  input 

function  edit4_Callback(h0bject ,  eventdata,  handles) 
handles. R  =  str2double (get (hObj ect ,  ’ String' )) ; 
guidata(hObject .handles) ; 

function  edit4_CreateFcn(h0bject ,  eventdata,  handles) 
handles. R  =  str2double (get (hObj ect ,  ’ String' )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObj ect , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j ect , ’ BackgroundColor ’ , J  white ’ ) ; 

end 

°/„Create  temperature  input 

function  edit5_Callback(h0bject ,  eventdata,  handles) 
handles. T  =  str2double (get (hObj ect , ’ String' )) ; 
guidata(hObject .handles) ; 

function  edit5_CreateFcn(h0bject ,  eventdata,  handles) 
handles. T  =  str2double (get (hObj ect String' )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObj ect , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

°/„Create  molecular  weight  input 

function  edit6_Callback(h0bject ,  eventdata,  handles) 
handles. M  =  str2double (get (hObj ect String' )) ; 
guidata(hObject .handles) ; 

function  edit6_CreateFcn(h0bject ,  eventdata,  handles) 
handles. M  =  str2double (get (hObj ect String' )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObj ect , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

7»Create  ambient  pressure  input 

function  edit7_Callback(h0bject ,  eventdata,  handles) 
handles. Po  =  str2double (get (hObj ect, ’String')) ; 
guidata(hObject .handles) ; 

function  edit7_CreateFcn(h0bject ,  eventdata,  handles) 
handles. Po  =  str2double (get (hObj ect, ’String')) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObj ect , ’BackgroundColor ’) , 
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get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

“/Create  kinematic  viscosity  input 

function  edit8_Callback(h0bject ,  eventdata,  handles) 
handles. v  =  str2double (get (hObj ect ,  ’ String' )) ; 
guidata(hObject .handles) ; 

function  edit8_CreateFcn(h0bject ,  eventdata,  handles) 
handles. v  =  str2double (get (hObj ect ,  ’ String' )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObj ect , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j ect , ’ BackgroundColor ’ , ’  white ’ ) ; 

end 

“/Create  bulk  viscosity  ratio  input 

function  edit9_Callback(h0bject ,  eventdata,  handles) 
handles. u_B  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

function  edit9_CreateFcn(h0bject ,  eventdata,  handles) 
handles. u_B  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

“/Create  Prandtl  number  input 

function  edit 10_Callback(h0bj ect ,  eventdata,  handles) 
handles. Pr  =  str2double (get (hObject, ’String')) ; 
guidata(hObject .handles) ; 

function  editlO_CreateFcn(hObject ,  eventdata,  handles) 
handles. Pr  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

%This  code  runs  when  button  is  pushed 

function  pushbuttonl_Callback (hObject ,  eventdata,  handles) 
clc ; 


SPL  =  handles. SPL; 
f  =  handles . f ; 


“/.Input  sound  pressure  level 
“/frequency 
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Po  =  handles. Po; 

M  =  handles. M; 
tempC  =  handles. T; 
R  =  handles. R; 

Y  =  handles. Y; 
v  =  handles. v; 
u_B  =  handles. u_B; 
Pr  =  handles. Pr; 


°/0ambient  pressure 

"/molecular  mass  of  air 

"/temp  Celsius 

"/ideal  gas  constant 

"/ratio  of  specific  heat  for  air 

"/kinematic  viscosity  coefficient 

"/bulk  viscosity  coefficient 

"/Prandtl  number 


Rs  =  R/M; 

c  =  331  +  0.6*tempC; 
tempK  =  273.15  +  tempC; 
rho_o  =  Po/ (Rs*tempK) ; 

Zo  =  c*rho_o; 
w  =  2*pi*f; 

P  =  sqrt (2)*((10~ (SPL/20) ) *20e-6) ; 
U  =  P/ (rho_o*c) ; 

Mu  =  U/c; 
k  =  w/c; 

V  =  (4/3) +u_B ; 


"/specific  gas  constant 

"/speed  of  sound 

"/temp  kelvin 

"/ambient  density 

"/characteristic  impedence 

"/radial  frequency 

"/Conversion  to  pressure  amplitude 

"/Conversion  to  particle  veloctiy 

"/acoustic  mach  number 

"/wave  number 

"/viscosity  number 


A  =  rho_o* (c~2) ; 

"/1st  derivative  of  sound  speed  with  respect  to  pressure 
dcdp  =  ( ( ( (P+Po) /Po) ~ ( (Y-l) / (2*Y) ) ) * (Y-l) *c) / (2*Y* (P+Po) ) ; 
"/2nd  derivative  of  sound  speed  with  respect  to  pressure 
d2cdp2  =  ( ( ( ( (P+Po) /Po) ~ ( (Y-l) / (2*Y)))*((Y-l)"2)*c)/ (4* (Y~2) 
*  (  (P+Po) ~2) ))-(((( (P+Po)/Po) ~ ( (Y-l) / (2*Y) ) ) 

* (Y-l) *c) / (2*Y* ( (P+Po) "2) ) ) ; 

BoA  =  2*rho_o*c*dcdp; 

CoA  =  ( (3/2) * (BoA~2) ) + (2* (rho_o~2) * (c~3) *d2cdp2) ; 

B  =  BoA*A; 

C  =  CoA*A; 


Bair  =  l+(0.5*(BoA)) ;  °/2nd  order  nonlinear  coefficient 

Nair  =  l+(0.5*(CoA)) ;  °/3rd  order  nonlinear  coefficient 

Gair  =  pi/ (2* (rho_o~2) * (c~5) ) ; 

Fair  =  6*pi/c; 

x_D  =  1/ (k*Bair*Mu) ;  "/discontinuity  distance 

"/1st  order  attenuation  coefficient 

alpha_l  =  ( ( (l*w) ~2) *v* (V+( (Y-l) /Pr) ) ) / (2* (c"3) ) ; 

"/2nd  order  attenuation  coefficient 

alpha_2  =  ( ( (2*w) ~2) *v* (V+( (Y-l) /Pr) ) ) / (2* (c"3) ) ; 
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0/03rd  order  attenuation  coefficient 

alpha_3  =  ( ( (3*w) ~2) *v* (V+( (Y-l) /Pr) ) ) / (2* (c~3) ) ; 


Fs  =  500000; 
time_samps  =  50; 
dt  =  1/Fs ; 

t  =  0 : dt : (time_samps-l) *dt ; 
x  =  x_D/ 1000 : x_D/ 1000 : 2*x_D ; 


%sample  rate 
%number  of  time  samples 
%sample  period 
%time  vector 

%1  dimension  propagation  distance 


for  i=l : length(x) 

"/(perturbation  equations  for  1st,  2nd,  and  3rd  harmonic 
P 1  ( i ,  : )  =  exp(-alpha_l*x(i))*P*sin(w*t-k*x(i)) ; 

P2 ( i , : )  =  exp(-alpha_2*x(i))*((Bair*w*(P~2) 

*x(i) ) / (2*rho_o* (c~3)))*sin(2* (w*t-k*x(i) ) ) ; 
P3(i,:)  =  exp(-alpha_3*x(i))*(Gair*f*(P"3)) 

* ( (Fair*f * (Bair" 2) * (x(i) ~2) )+(2*Nair*x(i) ) ) 

*sin(3* (w*t-k*x(i) ) ) ; 


°/„bessel  equations  for  1st,  2nd,  and  3rd  harmonic 
B 1  ( i , : )  =  2*P*(besselj (l,x(i)/x_D)/(x(i)/x_D)) 
*sin(l* (w*t-k*x(i) ) ) ; 

B2(i, : )  =  2*P*(besselj (2,2*x(i)/x_D)/(2*x(i)/x_D)) 
*sin(2* (w*t-k*x(i) ) ) ; 

B3(i,:)  =  2*P*(besselj (3,3*x(i)/x_D)/(3*x(i)/x_D)) 
*sin(3* (w*t-k*x(i) ) ) ; 


/(ideal  gas  treatment  for  1st,  2nd,  and  3rd  harmonic 
SI (i ,  : )  =  exp(-alpha_l*x(i))*P*sin(w*t-k*x(i)) ; 

S2 ( i , : )  =  exp(-alpha_2*x(i))*(P*0.5*Mu*Bair*k*x(i)) 

. *sin(2* (w*t-k*x(i) ) ) ; 

S3 ( i , : )  =  exp(-alpha_3*x(i))*(P*(Mu~2)*((Bair*k*x(i)) . ~2)) 
. * ( ( (sin(w*t-k*x(i) ) ) . * (cos (w*t-k*x(i) ) .  ~2) ) 
-(0.5*sin(w*t-k*x(i)) . ~3)) ; 

end 


°/0add  all  harmonic  sinusoids  together 
P_tot  =  P1+P2+P3; 

B_tot  =  B1+B2+B3; 

S_tot  =  S1+S2+S3 ; 

set (handles . textll , ’String’ ,alpha_l) ; 
set (handles. textl3, ’String’ ,alpha_2) ; 
set (handles . textl5 , ’String’ ,alpha_3) ; 
set (handles. textl7, ’String’ ,A) ; 


68 


set (handles . textl9 , ’String' ,B) ; 
set (handles . text21 , ’String’ ,C) ; 
set (handles . text29 , ’ String ’ , BoA) ; 
set (handles. text23, ’String’ ,  CoA) ; 
set (handles . text25 , ’ String ’ , x_D) ; 
set (handles . text31 , ’String’ ,Bair) ; 
set (handles. text33, ’String’ ,Nair) ; 

axes (handles . axesl) 
guidata(hObject .handles) ; 
subplot (3,1,1) 
hold  on 

plot (t ,P_tot (1 , : ) , ’k’ ) 
plot (t ,B_tot (1 , : ) , ’k: ’ ) 
plot (t ,S_tot (1 , : ) , ’k — ’ ) 
hold  off 

title (’Time  domain  waveform  for  x  =  O’); 
xlabel(’time  (seconds)’); 
ylabel (’ amplitude  (pascals)’); 

subplot (3,1,2) 
hold  on 

plot (t ,P_tot (1001 , :) , ’k’) 
plot (t,B_tot (1001, : ) , ’k: ’ ) 
plot (t , S_tot (1001 ,  :) ,  ’k— ’) 
hold  off 

legend ( ’Perturbation’ , ’Bessel-Fubini ’ , ’Ideal  Gas’) 
title(’Time  domain  waveform  for  x  =  x_D’); 
xlabel(’time  (seconds)’); 
ylabel (’ amplitude  (pascals)’); 

subplot (3,1,3) 
hold  on 

plot (t ,P_tot (2000 , : ) , ’k’ ) 
plot (t , B_tot (2000 , : ) , ’ k : ’ ) 
plot  (t ,  S_tot  (2000 ,  : ) ,  ’  k—  ’ ) 
hold  off 

title(’Time  domain  waveform  for  x  =  2x_D’); 
xlabel(’time  (seconds)’); 
ylabel (’ amplitude  (pascals)’); 
guidata(h0bject .handles) ; 

for  i=l : length(x) 

/(Perturbation  amplitudes 
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PI  (i , : )  =  P*exp(-alpha_l*x(i) ) ; 

P2(i, : )  =  ( ( (Bair*w* (P~2) *x(i) )/ (2*rho_o* (c~3)))) 
*exp(-alpha_2*x(i) ) ; 

P3(i,:)  =  ( (Gair*f * (P~3) ) * ( (Fair*f * (Bair~2) * (x(i) ~2) ) 
+(2*Nair*x(i) ) ) ) *exp(-alpha_3*x(i) ) ; 

%Bessel  amplitudes 

B 1  ( i ,  : )  =  2*P*(besselj (l,x(i)/x_D)/(x(i)/x_D)) ; 

B2  (i , : )  =  2*P*(besselj (2,2*x(i)/x_D)/(2*x(i)/x_D)) ; 

B3(i,:)  =  2*P*(besselj (3,3*x(i)/x_D)/(3*x(i)/x_D)) ; 

°/„Ideal  gas  amplitudes 
S 1  ( i ,  : )  =  P*exp(-alpha_l*x(i) ) ; 

S2(i,:)  =  ( (P*0 . 5*Mu*Bair*k*x(i) ) ) *exp(-alpha_2*x(i) ) ; 

S3(i , : )  =  (0.5* (P* (Mu~2) * ( (Bair*k*x(i) ) . ~2) ) ) *exp (-alpha_3*x (i) ) ; 

end 

Pl_dB  =  20*logl0 ( (Pl/sqrt (2) ) /20e-6) ; 

P2_dB  =  20*logl0((P2/sqrt(2))/20e-6) ; 

P3_dB  =  20*logl0((P3/sqrt(2))/20e-6) ; 

Bl_dB  =  20*logl0 ( (Bl/sqrt (2) ) /20e-6) ; 

B2_dB  =  20*logl0((B2/sqrt(2))/20e-6) ; 

B3_dB  =  20*logl0((B3/sqrt(2))/20e-6) ; 

Sl_dB  =  20*logl0 ( (Sl/sqrt (2) ) /20e-6) ; 

S2_dB  =  20*logl0((S2/sqrt(2))/20e-6) ; 

S3_dB  =  20*logl0((S3/sqrt(2))/20e-6) ; 

axes (handles . axes3) 
guidata(hObject ,  handles) 
subplot (1 , 3,1) 
hold  on 

plot (x , Pl_dB ,  ’  k  ’ ) ; 
plot(x,Bl_dB, ’k: ’) ; 
plot(x,Sl_dB, ’  k — ’) ; 
hold  off 

axis([0  2*x_D  0  130]) 

subplot (1,3, 2) 
hold  on 

plot (x , P2_dB ,  ’  k  ’ ) ; 
plot (x , B2_dB , ’ k : ’) ; 
plot (x , S2_dB , ’ k — ’ ) ; 
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hold  off 

axis([0  2*x_D  0  130]) 

legend ( ’Perturbation’ , ’Bessel-Fubini ’ , ’Ideal  Gas’) 

subplot (1,3, 3) 
hold  on 

plot(x,P3_dB, ’k’) ; 
plot (x , B3_dB , ’ k : ’) ; 
plot (x , S3_dB , ’ k — ’ ) ; 
hold  off 

axis([0  2*x_D  0  130]) 
guidata(hObject ,  handles) 

'/.sum  the  RMS  pressure  of  each  harmonic  together 

P_tot  =  (PI/ sqrt (2) ) .~2+ (P2/ sqrt (2) ) ,~2+ (P3/ sqrt (2) ) . ~2 ; 

B_tot  =  (Bl/sqrt(2)) . "2+ (B2/sqrt (2) ) . ~2+ (B3/sqrt (2) ) . ~2; 

S_tot  =  (SI/ sqrt (2)) . ~2+ (S2/ sqrt (2) ) . ~2+ (S3/ sqrt (2) ) . ~2 ; 

"/compute  total  pressure  in  dB  scale 

P_tot_dB  =  10*logl0(P_tot/(20e-6)~2) ; 

B_tot_dB  =  10*logl0 (B_tot/ (20e-6) ~2) ; 

S_tot_dB  =  10*logl0(S_tot/(20e-6)~2) ; 

axes (handles . axes4) 
guidata(hObject ,  handles) 
hold  on 

plot (x , P_tot_dB , ’ k ’ ) 
plot(x,B_tot_dB, ’k: ’) 
plot (x , S_tot_dB , ’ k — ’ ) 
hold  off 

legend ( ’Perturbation’ , ’Bessel-Fubini’ , ’Ideal  Gas’) 
axis ( [0  2*x_D  100  130]) 
guidata(h0bject ,  handles) 
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Appendix  D.  MUSIC  GUI  Application  Matlab  Code 


function  varargout  =  MUSIC_GUI_2(varargin) 

"/.  MUSIC_GUI_2  M-file  for  MUSIC_GUI_2 . f ig 

7.  MUSIC_GUI_2,  by  itself,  creates  a  new  MUSIC_GUI_2  or  raises  the 

7o  existing  singleton*. 

1 

7o  H  =  MUSIC_GUI_2  returns  the  handle  to  a  new  MUSIC_GUI_2  or  the 

7o  handle  to  the  existing  singleton*. 

1 

°l  MUSIC_GUI_2( ’CALLBACK ’ ,hObject , eventData, handles ,... )  calls  the 

%  local  function  named  CALLBACK  in  MUSIC_GUI_2 .M  with  the  given 

7o  input  arguments.  MUSIC_GUI_2( ’Property’ , ’Value’ ,... )  creates  a  new 

7o  MUSIC_GUI_2  or  raises  the  existing  singleton*.  Starting  from  the 

7o  left,  property  value  pairs  are  applied  to  the  GUI  before 

7o  MUSIC_GUI_2_0peningFunction  gets  called.  An  unrecognized  property 

7o  name  or  invalid  value  makes  property  application  stop.  All  inputs 

7o  are  passed  to  MUSIC_GUI_2_0peningFcn  via  varargin.  *See  GUI 

7o  Options  on  GUIDE'S  Tools  menu.  Choose  "GUI  allows  only  one 

7o  instance  to  run  (singleton) " . 
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°/o 

/  See  also:  GUIDE,  GUIDATA,  GUIHANDLES 

*/.  Edit  the  above  text  to  modify  the  response  to  help  MUSIC_GUI_2 

%  Last  Modified  by  GUIDE  v2.5  26-Feb-2008  10:42:34 

*/.  Begin  initialization  code  -  DO  NOT  EDIT 
gui_Singleton  =  1; 

gui_State  =  struct( ,gui_Name) ,  mfilename,  ... 

’ gui_Singleton’ ,  gui_Singleton,  . . . 

’ gui_OpeningFcn’ ,  @MUSIC_GUI_2_0peningFcn,  ... 

’ gui_0utputFcn’ ,  @MUSIC_GUI_2_0utputFcn,  . . . 

’ gui_LayoutFcn’ ,  [],... 

’ gui_Callback’ ,  [] ) ; 

if  nargin  &&  ischar (varargin{l}) 

gui_State .  gui_Callback  =  str2func  (varargin{l]-) ; 

end 

if  nargout 

[varargout{l :nargout}]  =  gui_mainf cn(gui_State ,  varargin{ : }) ; 

else 

gui_mainf cn(gui_State ,  varargin{ : }) ; 

end 

*/.  End  initialization  code  -  DO  NOT  EDIT 

l  -  Executes  just  before  MUSIC_GUI_2  is  made  visible. 

function  MUSIC_GUI_2_0peningFcn(h0bject ,  eventdata,  handles,  varargin) 
l  This  function  has  no  output  args,  see  OutputFcn. 

I  hObject  handle  to  figure 

l  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
l  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

*/.  varargin  command  line  arguments  to  MUSIC_GUI_2  (see  VARARGIN) 

*/.  Choose  default  command  line  output  for  MUSIC_GUI_2 
handles . output  =  hObject; 

*/.  Update  handles  structure 
guidata (hObject ,  handles); 

%  UIWAIT  makes  MUSIC_GUI_2  wait  for  user  response  (see  UIRESUME) 

*/.  uiwait  (handles .  figurel)  ; 

*/. - Outputs  from  this  function  are  returned  to  the  command  line. 
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function  varargout  =  MUSIC_GUI_2_0utputFcn(h0bject,  eventdata,  handles) 
%  varargout  cell  array  for  returning  output  args  (see  VARARGOUT) ; 

*/.  hObject  handle  to  figure 

*/.  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
*/.  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

*/.  Get  default  command  line  output  from  handles  structure 
varargout {1}  =  handles . output ; 

°/„Create  frequency  fl  input 

function  editl_Callback(hObject ,  eventdata,  handles) 
handles.fi  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

function  editl_CreateFcn(hObject ,  eventdata,  handles) 
handles.fi  =  str2double (get (hObject ,’ String’ )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

°/„Create  frequency  start  input 

function  edit2_Callback (hObject ,  eventdata,  handles) 

handles . freq_sweep_st  =  str2double (get (hObject ,’ String’ )) ; 
guidata(hObject .handles) ; 

function  edit2_CreateFcn(h0bject ,  eventdata,  handles) 

handles . freq_sweep_st  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

°/„Create  frequency  increment  input 

function  edit3_Callback (hObject ,  eventdata,  handles) 

handles . delta_f  =  str2double (get (hObject ,’ String’ )) ; 
guidata(hObject .handles) ; 

function  edit3_CreateFcn(h0bject ,  eventdata,  handles) 

handles . delta_f  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

°/„Create  frequency  sweep  stop  input 

function  edit4_Callback (hObject ,  eventdata,  handles) 

handles . freq_sweep_sp  =  str2double (get (hObject ,’ String’ )) ; 


guidata(hObject .handles) ; 

function  edit4_CreateFcn(h0bject ,  eventdata,  handles) 

handles . freq_sweep_sp  =  str2double (get (hObject, ’String')) 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  sampling  frequency  input 

function  edit5_Callback(h0bject ,  eventdata,  handles) 
handles. Fs  =  str2double (get (hObject, ’String')) ; 
guidata(hObject .handles) ; 

function  edit5_CreateFcn(h0bject ,  eventdata,  handles) 
handles. Fs  =  str2double (get (hObject ,’ String’ )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

7»Create  time  samples  input 

function  edit6_Callback (hObject ,  eventdata,  handles) 

handles . time_samps  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

function  edit6_CreateFcn(h0bject ,  eventdata,  handles) 

handles . time_samps  =  str2double (get (hObject ,’ String’ )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  input  file  name  input 

function  edit7_Callback (hObject ,  eventdata,  handles) 

handles . Hollow_Target_File  =  get (hObject ,’ String’ ) ; 
guidata(hObject .handles) ; 

function  edit7_CreateFcn(h0bject ,  eventdata,  handles) 
handles . Hollow_Target_File  =  get (hObject, ’String’) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

7»Create  input  file  name  input 

function  edit8_Callback (hObject ,  eventdata,  handles) 
handles . Solid_Target_File  =  get (hObject, ’String’ ) ; 


guidata(hObject .handles) ; 

function  edit8_CreateFcn(h0bject ,  eventdata,  handles) 
handles . Solid_Target_File  =  get (hObject, ’String’) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  output  file  name  input 

function  edit9_Callback(h0bject ,  eventdata,  handles) 
handles . Output_File  =  get (hObject ,  ’ String' ) ; 
guidata(hObject .handles) ; 

function  edit9_CreateFcn(h0bject ,  eventdata,  handles) 
handles . Output_File  =  get (hObject, ’String’) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  Fitlter  LSE  order  input 

function  editl8_Callback(h0bject ,  eventdata,  handles) 

handles . LSE_0rder  =  str2double (get (hObject, ’String')) ; 
guidata(hObject .handles) ; 

function  editl8_CreateFcn(h0bject ,  eventdata,  handles) 
handles . LSE_0rder  =  str2double (get (hObject, ’String')) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  Filter  time  samples  input 

function  editlO_Callback(hObject ,  eventdata,  handles) 

handles . Filter_Samples  =  str2double (get (hObject ,’ String’ )) ; 
guidata(hObject .handles) ; 

function  editlO_CreateFcn(hObject ,  eventdata,  handles) 

handles . Filter_Samples  =  str2double (get (hObject, ’String’)) 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

7»Create  Filter  start  frequency  input 

function  editll_Callback(hObject ,  eventdata,  handles) 

handles . Filter_Start_Freq  =  str2double (get (hObject, ’String’)) ; 
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guidata(hObject .handles) ; 

function  editll_CreateFcn(hObject ,  eventdata,  handles) 

handles . Filter_Start_Freq  =  str2double (get (hObject , ’ String’ )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/.Create  Filter  stop  frequency  input 

function  editl2_Callback(h0bject ,  eventdata,  handles) 

handles . Filter_Stop_Freq  =  str2double (get (hObject ,  ’ String’ )) ; 
guidata(hObject .handles) ; 

function  editl2_CreateFcn(h0bject ,  eventdata,  handles) 

handles . Filter_Stop_Freq  =  str2double (get (hObject , ’ String’ )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/.Create  Filter  plot  maximum  y-axis  input 

function  edit20_Callback (hObject ,  eventdata,  handles) 

handles . Filter_Max  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

function  edit20_CreateFcn (hObject ,  eventdata,  handles) 
handles . Filter_Max  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/.Create  Filter  plot  minimum  y-axis  input 

function  edit21_Callback(h0bject ,  eventdata,  handles) 

handles . Filter_Min  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

function  edit21_CreateFcn(h0bject ,  eventdata,  handles) 
handles . Filter_Min  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/.Create  MUSIC  plot  maximum  y-axis  input 

function  edit22_Callback (hObject ,  eventdata,  handles) 

handles .Music_Max  =  str2double (get (hObject, ’String’)) ; 


guidata(hObject .handles) ; 

function  edit22_CreateFcn(h0bject ,  eventdata,  handles) 
handles .Music_Max  =  str2double(get(h0bject, ’String')) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’ ) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  MUSIC  plot  minimum  y-axis  input 

function  edit23_Callback (hObject ,  eventdata,  handles) 

handles .Music_Min  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

function  edit23_CreateFcn(h0bject ,  eventdata,  handles) 
handles .Music_Min  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

7»Create  MUSIC  time  samples  input 

function  edit24_Callback (hObject ,  eventdata,  handles) 

handles .music_samps  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

function  edit24_CreateFcn (hObject ,  eventdata,  handles) 
handles .music_samps  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/oCreate  MUSIC  number  of  sinusoids  to  detect  input 
function  edit 14_Callback(h0bj ect ,  eventdata,  handles) 

handles .num_of_sins  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

function  editl4_CreateFcn(h0bject ,  eventdata,  handles) 

handles .num_of_sins  =  str2double (get (hObject ,’ String’ ) ) 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

7»Create  MUSIC  segment  length  input 

function  editl5_Callback(h0bject ,  eventdata,  handles) 

handles . segment_length  =  str2double (get (hObject ,’ String’ )) ; 
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guidata(hObject .handles) ; 

function  editl5_CreateFcn(h0bject ,  eventdata,  handles) 

handles . segment_length  =  str2double (get (hObject , ’String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/.Create  MUSIC  overlap  percent  input 

function  editl6_Callback(h0bject ,  eventdata,  handles) 

handles . overlap_percent  =  str2double (get (hObject , ’ String’ )) ; 
guidata(hObject .handles) ; 

function  editl6_CreateFcn(h0bject ,  eventdata,  handles) 

handles . overlap_percent  =  str2double (get(hObj ect, 'String’)) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/.Create  MUSIC  threshold  input 

function  editl7_Callback(h0bject ,  eventdata,  handles) 

handles . threshold  =  str2double (get (hObject, ’String’)) ; 
guidata(hObject .handles) ; 

function  editl7_CreateFcn(h0bject ,  eventdata,  handles) 

handles . threshold  =  str2double (get (hObject ,’ String’ )) ; 
guidata(hObject .handles) ; 

if  ispc  &&  isequal (get (hObject , ’BackgroundColor ’) , 
get (0, ’ def aultUicontrolBackgroundColor ’ ) ) 
set (hOb j  ect , ’ BackgroundColor ’ , ’ white ’ ) ; 

end 

"/.Excecute  on  button  press 

function  pushbuttonl_Callback (hObject ,  eventdata,  handles) 
clc 

guidata(hObject .handles) ; 
try 

°/„read  inputs  from  GUI 
f req_sweep_st  =  handles . freq_sweep_st ; 
delta_f  =  handles . delta_f ; 
f req_sweep_sp  =  handles . freq_sweep_sp; 

f2  =  freq_sweep_st:delta_f :freq_sweep_sp; 

fl  =  linspace (handles . f 1 .handles . fl , length (f 2) ) ; 

d_f  =  f 2-f 1 ; 


Fs  =  handles. Fs; 

time_samps  =  handles. time_samps; 
dt  =  1/Fs; 

Order  =  handles. LSE_Order; 
filt_samps  =  handles . Filter_Samples ; 
filt_st  =  handles . Filter_Start_Freq; 
filt_sp  =  handles . Filter_Stop_Freq; 

music_samps  =  handles. music_samps; 

NS  =  handles .num_of_sins ; 

SL  =  handles . segment_length; 

OLP  =  handles . overlap_percent ; 

Thresh  =  handles .threshold; 


°/0load  data 

Hollow_Data  =  load(char (handles .Hollow_Target_File) ) ; 
Solid_Data  =  load(char (handles . Solid_Target_File) ) ; 
"/parse  data 
try 

for  ind=l : length (d_f) 

l_ind  =  ((ind-l)*time_samps)+l; 
u_ind  =  ind*time_samps ; 

Hollow_inc(ind, : )  =  Hollow_Data(l_ind:u_ind,2) 
Hollow_ref (ind, : )  =  Hollow_Data(l_ind:u_ind,3) 
Solid_inc(ind, : )  =  Solid_Data(l_ind:u_ind,2) ; 
Solid_ref (ind, : )  =  Solid_Data(l_ind:u_ind,3) ; 


end 

"/initialize  output  values 
Z_HI  =  0; 

Z_HR  =  0; 

Z_SI  =  0; 

Z_SR  =  0; 


M_f req_HI 
M_f req_HR 
M_f req_SI 
M_f req_SR 


zer os (1, length (d_f)) ; 
zer os (1, length (d_f)) ; 
zer os (1, length (d_f)) ; 
zer os (1, length (d_f)) ; 


M_amp_HI  =  zeros (1 , length(d_f )) 
M_amp_HR  =  zeros (1 , length(d_f )) 
M_amp_SI  =  zeros (1 , length(d_f )) 
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M_amp_SR  =  zeros (1 , length(d_f )) ; 


M_amp_HI_dB  =  zeros(l , length(d_f ) ) ; 
M_amp_HR_dB  =  zeros(l , length(d_f ) ) ; 
M_amp_SI_dB  =  zeros(l , length(d_f ) ) ; 
M_amp_SR_dB  =  zeros (1 , length (d_f) ) ; 


for  ind=l : length (d_f) 

FOI  =  ind*delta_f ;  °/„Frequency  of  interest 
set (handles. text 44, 'String’ ,ind) ; 
set (handles. text 43, ’String’ ,F0I) ; 

7»read  block  of  noisy  data 

S_n_t  =  Hollow_inc(ind, 1 : f ilt_samps) ; 

fft_st  =  0; 

fft_sp  =  d_f (length(d_f ) ) ; 

°/oCall  the  prewhitening  filter 
[S_w_t,  S_w_f ,  S_n_f ,  freq,  apx_dB, 
filt_freq,  CC]  =  f_white_GUI 
(S_n_t(l :filt_samps) ,dt,fft_st,fft_sp, 
filt_st,filt_sp,ind,delta_f ,  Order) ; 

7»plot  noisy  data,  approximation,  and  whitening 
7»data 

axes (handles . axes3) 
cla (handles . axes3 , ’reset ’ ) ; 
guidata(hObject, handles) ; 
hold  on 

plot (f req, S_n_f , ’r ’ ) ; 
plot (f req, S_w_f , ’b’ ) ; 
plot (filt_f req,  apx_dB,  ’g’); 
axis ( [f f t_st , f f t_sp , handles . Filter_Min, 
handles .Filter_Max] ) ; 
hold  off 

set (handles . text20 , ’ String ’ , CC) ; 
guidata(hObject, handles) ; 

/(analyze  data  with  MUSIC  algorithm 
d  =  (floor ( (f ilt_samps-music_samps) /2) )+l ; 
[M_amp_HI (ind) ,  M_f req_HI (ind) ,  Spectrum, 
Frequencies,  Z]  =  my_music_GUI 
(S_w_t(d: (d-l)+music_samps) ,Fs,F0I, 
delta.f,  NS,  SL,  OLP,  Thresh); 

7»keep  track  of  zeroed  coefficients 
°/0and  plot  data 
Z_HI  =  Z_HI  +  Z; 

P_HI  =  100* (1- (Z_HI/length(d_f ) ) ) ; 
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axes (handles . axes4) 
cla (handles . axes4, ’reset ’ ) ; 
guidata(hObject, handles) ; 
hold  on 

plot (Frequencies , Spectrum, ’b’ ) ; 
stem(M_f req_HI (ind) ,M_amp_HI (ind) , ’r ’ ) ; 
hold  off 

axis([0  10000  handles .Music_Min 
handles .Musi c_Max] ) ; 
set (handles. textl9, ’String’ ,Z_HI) ; 
set (handles. textl8, ’String’ ,P_HI) ; 
guidata(hObject, handles) ; 

S_n_t  =  Hollow_ref (ind, 1 : f ilt_samps) ; 
fft_st  =  0; 

fft_sp  =  d_f (length (d_f) ) ; 

[S_w_t,  S_w_f ,  S_n_f ,  freq,  apx_dB,  filt_freq, 
CC]  =  f _white_GUI (S_n_t (1 : f ilt_samps) ,dt , 
fft_st,fft_sp,filt_st , f ilt_sp, ind, delta_f , 
Order) ; 

axes (handles . axes5) 
cla (handles . axes5 , ’reset ’ ) ; 
guidata(h0bject, handles) ; 
hold  on 

plot (f req, S_n_f , ’r ’ ) ; 
plot (f req, S_w_f , ’b’ ) ; 
plot (filt_f req,  apx_dB,  ’g’); 
axis ( [f f t_st , f f t_sp, handles . Filter_Min, 
handles .Filter_Max] ) ; 
hold  off 

set (handles . text27 , ’ String ’ , CC) ; 
guidata(h0bject, handles) ; 

d  =  (floor ( (f ilt_samps-music_samps) /2) )+l ; 
[M_amp_HR(ind) ,  M_f req_HR(ind) ,  Spectrum, 
Frequencies,  Z]  =  my_music_GUI 
(S_w_t(d: (d-l)+music_samps) ,Fs,F0I, 
delta.f,  NS,  SL,  0LP,  Thresh); 

Z_HR  =  Z_HR  +  Z; 

P_HR  =  100* (1- (Z_HR/length(d_f ) ) ) ; 
axes (handles . axes6) 
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cla (handles . axes6 , ’reset ’ ) ; 
guidata(hObject,handles) ; 
hold  on 

plot (Frequencies , Spectrum, ’b’ ) ; 
stem(M_f req_HR(ind) ,M_amp_HR(ind) , ’r ’ ) ; 
hold  off 

axis([0  10000  handles .Music_Min 
handles .Musi c_Max] ) ; 
set (handles. text26, ’String’ ,Z_HR) ; 
set (handles. text25, ’String’ ,P_HR) ; 
guidata(hObject, handles) ; 

S_n_t  =  Solid_inc(ind, 1 :f ilt_samps) ; 
fft_st  =  0; 

fft_sp  =  d_f (length(d_f ) ) ; 

[S_w_t,  S_w_f ,  S_n_f ,  freq,  apx_dB,  filt_freq, 
CC]  =  f _white_GUI (S_n_t (1 : f ilt_samps) ,dt , 
fft_st,fft_sp,filt_st , f ilt_sp, ind, 
delta_f,  Order); 

axes (handles . axes7) 
cla (handles . axes7 , ’reset ’ ) ; 
guidata(h0bject, handles) ; 
hold  on 

plot (f req, S_n_f , ’r ’ ) ; 
plot (f req, S_w_f , ’b’ ) ; 
plot (filt_f req,  apx_dB,  ’g’); 
axis ( [f f t_st , f f t_sp, handles . Filter_Min, 
handles .Filter_Max] ) ; 
hold  off 

set (handles . text33 , ’ String ’ , CC) ; 
guidata(h0bject, handles) ; 

d  =  (floor ( (f ilt_samps-music_samps) /2) )+l ; 
[M_amp_SI (ind) ,  M_f req_SI (ind) ,  Spectrum, 
Frequencies,  Z]  =  my_music_GUI 
(S_w_t(d: (d-l)+music_samps) ,Fs,F0I, 
delta.f,  NS,  SL,  0LP,  Thresh); 

Z_SI  =  Z_SI  +  Z; 

P_SI  =  100* (1- (Z_SI/length(d_f ) ) ) ; 

axes (handles . axes8) 

cla (handles . axes8 , ’reset ’ ) ; 
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guidata(hObject,handles) ; 
hold  on 

plot (Frequencies , Spectrum, ’b’ ) ; 
stem(M_f req_SI (ind) ,M_amp_SI (ind) , ’r ’ ) ; 
hold  off 

axis([0  10000  handles .Music_Min 
handles .Musi c_Max] ) ; 
set (handles. text32, 'String’ ,Z_SI) ; 
set (handles. text31, ’String’ ,P_SI) ; 
guidata(hObject, handles) ; 

S_n_t  =  Solid_ref (ind, 1 :f ilt_samps) ; 
fft_st  =  0; 

fft_sp  =  d_f (length(d_f ) ) ; 

[S_w_t,  S_w_f ,  S_n_f ,  freq,  apx_dB,  filt_freq, 
CC]  =  f _white_GUI (S_n_t (1 : f ilt_samps) ,dt , 
fft_st,fft_sp,filt_st , f ilt_sp, ind, 
delta_f,  Order); 

axes (handles . axes9) 
cla (handles . axes9 , ’reset ’ ) ; 
guidata(h0bject, handles) ; 
hold  on 

plot (f req, S_n_f , ’r ’ ) ; 
plot (f req, S_w_f , ’b’ ) ; 
plot (f ilt_freq,  apx_dB,  ’g’); 
axis ( [f f t_st , fft_sp, handles . Filter_Min, 
handles .Filter_Max] ) ; 
hold  off 

set (handles . text39 , ’ String ’ , CC) ; 
guidata(h0bject, handles) ; 

d  =  (floor ( (f ilt_samps-music_samps) /2) )+l ; 
[M_amp_SR(ind) ,  M_f req_SR(ind) ,  Spectrum, 
Frequencies,  Z]  =  my_music_GUI 
(S_w_t(d: (d-l)+music_samps) ,Fs,F0I, 
delta.f,  NS,  SL,  0LP,  Thresh); 

Z_SR  =  Z_SR  +  Z; 

P_SR  =  100* (1- (Z_SR/length(d_f ) ) ) ; 
axes (handles . axes  10) 
cla (handles . axes 10 , ’reset ’ ) ; 
guidata(h0bject, handles) ; 
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hold  on 

plot (Frequencies , Spectrum, ’b’ ) ; 
stem(M_f req_SR(ind) ,M_amp_SR(ind) , ’r ’ ) ; 
hold  off 

axis([0  10000  handles .Music_Min 
handles .Musi c_Max] ) ; 
set (handles. text38, ’String’ ,Z_SR) ; 
set (handles. text37, ’String’ ,P_SR) ; 
guidata(hObject, handles) ; 

/.Convert  MUSIC  coefficients  to  dB  scale 
M_amp_HI_dB (ind)  =  20*logl0((M_amp_HI(ind)/ 
sqrt (2) ) /20e-6) ; 

M_amp_HR_dB (ind)  =  20*logl0((M_amp_HR(ind)/ 
sqrt (2) ) /20e-6) ; 

M_amp_SI_dB (ind)  =  20*logl0((M_amp_SI(ind)/ 
sqrt (2) ) /20e-6) ; 

M_amp_SR_dB (ind)  =  20*logl0((M_amp_SR(ind)/ 
sqrt (2) ) /20e-6) ; 

"/.plot  in  real  time 
axes (handles . axes  11) 
cla(handles . axesll , ’reset ’ ) ; 
guidata(hObject, handles) ; 
hold  on 

plot (M_f req_HI (1 : ind) ,M_amp_HI_dB(l : ind) , 
’b’); 

plot (M_freq_HR(l : ind) ,M_amp_HR_dB(l : ind) , 
’r’); 

plot (M_freq_SI (1 : ind) ,M_amp_SI_dB(l : ind) , 
’b: ’); 

plot (M_freq_SR(l : ind) ,M_amp_SR_dB(l : ind) , 
’r: ’); 
hold  off 

legend( ’Hollow  Incident’, 

’Hollow  Ref lected’ ,’ Solid  Incident’, 
’Solid  Reflected’) ; 
guidata(hObject .handles) ; 

end 

'/.save  output  data 

Output  =  [M_freq_HI’  M_amp_HI ’  M_freq_HR’  M_amp_HR’ 
M_f req_SI ’  M_amp_SI’  M_freq_SR’  M_amp_SR’] ; 
File_I0  =  handles . Output _File ; 
save  File_I0  Output; 

/throw  error  if  data  is  not  parsed  correctly 
catch 
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disp(’ Error  parsing  data  -  Please  check  frequency 
increments/time  samples’); 
lasterror 

end 

°/0throw  error  if  filenames  are  misspelled 
catch 

disp(’Error  loading  file(s)  -  Please  check  spelling/syntax’ ) ; 
lasterror 

end 

"/throw  error  if  non-numeric  data  is  entered 
catch 

disp( ’Error  reading  data  -  Please  check  frequency 
increments/time  samples’); 
lasterror 

end 


86 


Appendix  E.  Prewhitening  Filter  Matlab  Code 

function  [S_w_C , amp_f ilt_dB_C , amp_dB , f req, apx_dB , f ilt_f req, CC] 

=  f _white_GUI (S_n_t ,dt , f f t_st ,f f t_sp,f ilt_st , f ilt_sp, ind,delta_f ,  order) 


*/.  Prewhitening  Filter:  whitens  narrow  band  signals  in  1/f  noise 

% 


1 

1 

7. 

7. 

7 

7. 

*/. 

7. 

7 

7. 

7. 

7. 

7. 

% 

7. 

7o 

7. 

7. 

7o 

7o 


usage : 

f _white_GUI (S_n_t , dt , f f t_st , f f t_sp , f ilt_st , f ilt_sp , ind , delta_f , order) 

inputs:  S_n_t  =  noisy  time  domain  data 

dt  =  delta  t 

fft_st  =  DFT  start  frequency 
fft_sp  =  DFT  end  frequency 
filt_st  =  filter  start  frequency 
filt_sp  =  filter  stop  frequency 
ind  =  calling  loop  index 
delta_f  =  frequency  step 
order  =  LSE  approximation  order 

outputs:  S_w_C  =  whitened  time  domain  data 

amp_f ilt_dB_C  =  whitened  and  amplitude  corrected  spectrum 

amp_dB  =  noisy  unfiltered  spectrum 

freq  =  frequency  vector 

apx_dB  =  noise  approximation  in  dB 

filt_freq  =  filter  frequency  vector 

CC  =  correction  coefficient 


7«Take  the  DFT  of  the  input  signal 

[amp,freq,k,re, im]  =  mex_fft(S_n_t,dt,fft_st,fft_sp) ; 
7»Determine  the  frequency  resolution 
f_res  =  l/(dt*length(S_n_t) ) ; 

7oCompute  the  phase  component 
phi  =  atan(im/re) ; 

7oUse  dB  scale 

amp_dB  =  20*logl0((amp*(sqrt(2)/2))/20e-6) ; 
warning  ('off’); 

7oCompute  the  upper  and  lower  filter  frequencies 
l_ind  =  (f loor (f ilt_st/f _res) )+l ; 
u_ind  =  ceil(f ilt_sp/f_res) ; 
filt_freq  =  freq(l_ind:u_ind) ; 

7»Least  square  polyfit  of  noise  over  filter  frequencies 
P  =  polyfit (f ilt_freq,amp_dB(l_ind:u_ind) , order) ; 
warning  ( ’ on ’ ) ; 


"/.Approximation  of  noise  in  dB  and  its  minimum 
apx_dB  =  polyval(P,f ilt_freq) ; 
apx_min  =  min(apx_dB) ; 

"/.Some  code  to  deal  with  noise  floor  being  negative  dB 
if  apx_min<0 

apx_dB  =  -apx_dB; 
if  min(apx_dB) <0 

apx_dB  =  (apx_dB- (min(apx_dB) ) )+l ; 

end 

°/0The  filter  is  essentially  the  inverse  of  the  approximation 
filt_dB  =  l./apx_dB; 

°/0And  needs  to  be  normalized  to  one 
filt_dB  =  filt_dB*(l/max(filt_dB)) ; 
amp_filt_dB  =  l*amp_dB; 

amp_f ilt_dB(l : length (filt_dB) )  =  amp_dB(l : length(f ilt_dB) ) . *filt_dB 

else 

filt_dB  =  l./apx_dB; 

filt_dB  =  filt_dB*(l/max(filt_dB)) ; 

amp_filt_dB  =  l*amp_dB; 

amp_f ilt_dB(l : length (filt_dB) )  =  amp_dB(l : length(f ilt_dB) ) . *filt_dB 

end 

"/.Determine  what  the  correction  coefficient  is  and  what  index 
"/.it  applies  to 
if  (ind*delta_f )>=f ilt_sp 

CC_ind  =  ( ( (ind*delta_f )-filt_st)/f_res)+l; 

CC  =  1; 

else 

CC_ind  =  ( ( (ind*delta_f )-filt_st)/f_res)+l; 

CC  =  1/f ilt_dB (CC_ind) ; 

end 

freq(CC_ind) ; 

CC; 

"/.Rescale  the  filtered  amplitude  to  Pascals 

amp_filt  =  ( (10 . ~ (amp_f ilt_dB/20) ) *20e-6) / (sqrt (2) /2) ; 

"/.Find  the  real  and  imaginary  parts 
re_filt  =  amp_f ilt*cos (phi) ; 
im_filt  =  amp_f ilt*sin(phi) ; 

"/.Use  inverse  DFT  to  return  to  time  domain 

S_w  =  mex_ifft(re_filt,im_filt,k,length(S_n_t)) ; 

"/.Use  the  correction  coefficient 

amp_f ilt_dB_C  =  rescale (amp_f ilt_dB,  CC,  CC_ind) ; 


°/„Rescale  the  filtered  amplitude  to  Pascals 

amp_filt_C  =  ((10. "(amp_filt_dB_C/20))*20e-6)/(sqrt(2)/2) 

°/0Find  the  real  and  imaginary  parts 

re_filt_C  =  amp_f ilt_C*cos (phi) ; 

im_f ilt_C  =  amp_f ilt_C*sin(phi) ; 

7»Use  inverse  DFT  to  return  to  time  domain 

S_w_C  =  mex_if ft (re_f ilt_C , im_f ilt_C,k,length(S_n_t)) ; 
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Appendix  F.  MUSIC  Algorithm  Matlab  Code 

function  [M_amp,  M_freq,  M_spectrum,  M_f requencies ,  Z] 

=  my_music_GUI (S_w_t ,  Fs,  freq,  delta_f ,  NS,  SL,  OLP,  THRESH) 


*/.  Multiple  Signal  Classification 

1 


1  usage : 

°/o 

%  inputs: 

/ 

l 

l 

l 

1 

l 

1 

l 

°/„  outputs: 
•/. 

•/. 

i 

i 

i 

i 


my_music_GUI (S_w_t ,  Fs,  freq,  delta_f,  NS,  SL,  OLP  THRESH) 

S_w_t  =  whitened  time  domain  data 
Fs  =  sampling  frequency 

freq  =  the  current  frequency  being  analyzed 
delta_f  =  the  frequency  step  in  the  sweep 
NS  =  number  of  sinusoids  to  estimate 
SL  =  segment  length  used  to  partition  data 
OLP  =  overlap  percent  in  segmentation 

THRESH  =  the  threshold  between  noise  and  signal  subspace 

M_amp  =  the  amplitude  of  the  MUSIC  spectrum  at  the 
frequency  being  analyzed 

M_freq  =  the  frequency  estimate  of  the  MUSIC  spectrum  at 
the  frequency  being  analyzed 
M_spectrum  =  the  entire  MUSIC  spectrum  amplitudes 
M_f requencies  =  the  entire  MUSIC  spectrum  frequencies 
Z  =  flag  thrown  if  frequency  estimate  doesn’t  match 


°/0input  parameters 
WINNAME  =  ’Rectangular’ ; 

WINPARAM  =  ”  ; 

FFTLENGTH  =  ’NextPow2’; 

INPUTTYPE  =  ’Vector’; 

°/0define  the  pseudospectrum  estimator 

Hs  =  spectrum. music(NS,SL, OLP, WINNAME, THRESH, FFTLENGTH, INPUTTYPE) ; 


’/.use  Hs  to  calculate  the  pseudospectrum  of  S_w_t 

HPS  =  pseudospectrum (Hs , S_w_t , ’Fs ’ ,Fs , ’ SpectrumRange ’ , ’half ’ ) ; 

°/0Pull  out  the  amplitude  of  HPS 

M_spectrum  =  HPS. data; 

°/0Convert  to  dB 

M_spectrum  =  20*logl0((M_spectrum*sqrt(2)/2)/20e-6) ; 

°/0Pull  out  the  frequencies  of  HPS 
M_f requencies  =  HPS. freq; 


°/0Find  the  maximum  of  the  spectrum 
[M_amp,  ind]  =  max(M_spectrum) ; 
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M_freq  =  M_f requencies (ind) ; 

°/0Define  upper  and  lower  frequency  thresholds 
thresh_l  =  freq  -  (delta_f/2) ; 
thresh_u  =  freq  +  (delta_f/2) ; 

°/0See  if  MUSIC  frequency  estimate  is  within  threshold 
if  (thresh_l<M_freq  &&  M_f req<thresh_u) 

M_freq  =  l*M_freq; 

M_amp  =  l*M_amp; 

Z  =  0; 

7»if  it’s  not,  throw  the  flag, 
else 

M_freq  =  freq; 

°/0M_amp  gets  set  to  MUSIC  noise  floor  (30) 

M_amp  =  30; 

Z  =  1; 

end 
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Appendix  G.  Amplitude  Correction  Matlab  Code 

function  [amp_C]  =  rescale (amp,  CC,  CC_ind) 

%  Amplitude  Correction:  rescales  a  single  sinusoid  in  white  noise  to 
*/.  original  amplitude 
1 

l  usage:  rescale (amp,  CC,  CC_ind) 

l 

%  inputs:  amp  =  a  row  vector  containing  the  DFT  coefficients 

*/.  CC  =  the  computed  correction  coefficient  from 

%  ’f_white_GUI.m’ 

°/„  CC_ind  =  the  index  of  the  vector  ’amp’  that  needs  rescaling 

% 

l  outputs:  amp_C  =  the  rescaled  vector  of  DFT  coefficients 
amp_C  =  amp; 

amp_C(CC_ind)  =  CC*amp_C(CC_ind) ; 
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Appendix  H.  Equipment  Specifications 

National  Instruments  PXI-5922  High  Speed  Digitizer 

24-Bit  Flexible  Resolution  Digitizer 


Specifications 

These  specifications  a  re  valid  for  0°  to  55°  C.  unless  otherwise  stated.  Spectral  Characteristics  (typical) 


Acquisition  System 

Number  of  channels. . . . . . ...  2  simultaneously  sampled 

single-ended  or  unbalanced  differential 
or 

1  differential  channel 

Vertical  resolution _ _  1 G  to  24  bits 

Sampling  Rate  Resolutioi  I  bits) 

50  to  500  kS/s  24 

1  MS/s  22 

5  MS/s  20 

10  MS/s  18 

15  MS/s_ 18 


Alias-free  bandwidth  ..... _ _ _ _ _  0.4  x  sampling  rate 

Onboard  sample  memory - ....  8, 32.  or  25G  MB  per  channel  (2. 8.  or  64  million  samples) 

Pre  and  post  trigger  data  points- _  0-100%  of  full  record  length 

Multiple  Record  Acquisition  <0-100%  pre  aid  post  trigger  data) 

Memoiy  per  Channel  Maximum  Nundier  of  Records 

8  MB  13.107 

32  MB  52.428 

256  MB_ 419,430 


Input  impedance . . . . . .  50  O  and  1  Mil  II 40  pf  software  selectable 

Full  scale  input  range -  2  V^  |±1  V)  and  10  Vff  (±5  V) 

Maximum  input  overbad _  50  £17  with  |  peaks |  <  10  V  I  M£1  |  peaks |  <  42  V 

Input  coupling _ _ _ _ . _ .....  AC.  0C.  GNO 

AC  coupling  cutoff  frequency  (-3  dB) - -  90  Hz 

Accuracy 
DC  accuracy 

Range  50  £1  aid  1  MO 

2  Vpp  (±1  V)  ±(500  ppm  (0.05%)  of  Input  ±  50  pV) 

10  VPD  (±5  V) ±(500  ppm  (0.05%)  of  Input  1 100  pV) 


Note:  Measured  with  1  Mfl  input  impedance  within  ±5  °C  of  self-calibration  temperature. 

Passband  Flatness  (referenced  at  DC) 


Dynamic  performance  (-1  dBFS  input  signal) 


liput  Frequency 

SFDR(dBc) 

10VPp(±5VllnputRaiae 

THD(dBc) 

SINADldB) 

SNR  IrB) 

100  kHz 

-108 

-10G 

103 

106 

1  MHz 

-96 

-94 

89 

91 

liput  Frequency 

SFDR(dBc) 

2  Vpp  (±1  V)  Input  Range 
THD(dBc) 

SINADldB) 

SNR  l<B) 

100  kHz 

-103 

-101 

99 

103 

1  MHz 

-92 

-90 

87 

90 

Note:  Sampling  rate  is  1  MS/s  for  100  kHz  input  and  10  MS/s  for  1  MHz  input. 


Rms  Noise 

10  Vpp  (±5  V)  liput  Range  2  Vpp  (±1 V)  Inpit  Range 


Sample  Rate 

(dBFS) 

fllVrn.) 

(dBFS) 

50  kS/s 

-120 

3.4 

-117 

1.0 

lOOkS/s 

-118 

4.3 

-115 

1.2 

1  MS/s 

-108 

13 

-104 

4.2 

5  MS/s 

-101 

31 

-98 

8.7 

10  MS/s 

-91 

92 

-91 

20 

15  MS/s 

-79 

401 

-79 

80 

Phase  noise  density  |5  MHz  input' . - 

....  <-133  dBc/Hz  at  10  kHz. 

<-145  dBc/Hz  at  100  kHz 

Tune  base  System 

Timebase  options _  Internal 

Total  sample  clock  jitter _  ^3  ps^j 

Note:  Includes  effects  of  converter  aperture  and  clock  circuitry  jitter  from  100  Hz  to  1  MHz. 

Internal 

Internal  sample  clock  frequency -  1 20  MS/s  sampling  rate  with  decimation 

by  n  where  8  <  n  <  2400 

Timebase  accuracy,  typical _ _ . _ _  ±50  ppm  (±0.0050%) 

External 

External  reference  clock  sources . .  CLK  IN  (SMB  connector).  PXI  backplane  10  MHz 

External  reference  clock  range . . .  1  to  20  MHz  in  1  MHz  increments 


Sampling  Rate 

1  MS/s 

50  Cl  and  1  M£i,  ±1  V  and  ±5  V  ranges 
0.03  dB 

5  MS/s 

0.06  dB 

10  MS/s 

0.15  dB 

15  MS/s 

0.3  dB 

AC  amplitude  accuracy  (typical  at  1  kHz) .  ±  600  ppm  (0.06%) 

Channel  to  Channel  Crosstalk 

Input  Freqiency 

Crosstalc 

100  kHz 

< -110  dB 

1MHz 

< -100  dB 

G  MHz 

i-80  dB 

CMRR 

Unbalanced  differential  mode _ _  50  dB 


Differential  mode ..... _ .....  See  Figure  4 


Trigger  System 

Modes.. .  Edge,  hysteresis,  window,  digital,  immediate,  software 

Sources _ _ _ _ _ _ _  CH  0.  CH 1 .  TRIG.  PX)_Trig  <0f>>.  PFI  <0:1  >.  PXI  Star.  Software 

Slope -  Rising  or  falling 

Hysteresis- .  Fully  programmable 

Sensitivity _ ., _ , _ , _  CH  0  and  CH  1 : 2%  FS 

TRIG:  0.3  Vpp  typical  up  to  1MHz 

Time  resolution _ _ _ - . . .  One  sample  clock  period 

Rearm  time. _ _ _  144  x  sample  clock  period 

Hold  off  _ _ _ _ _ _ _ _ _ _ _  Up  to  |2K  -1)  x  sample  clock  period 

External  Trigger  Channel  (TRIG) 

Impedance _ - _ - _ _ _ _ _  100  Hi  II  52  pF 

Range  _  ±2.5  V 

Coupling..... . . . . . . . . . . . .  DC 

Level  Accuracy _ _ _  ±0.3  V  up  to  100  kHz 

Power  Requirements 

Voltage  (VDC)  Typical  Cuneat  IA) 

±3.3  2 

±5  1.4 

±12  0.33 

-12  028 

Total  Power_ 20.9  W 

Calibration 

Self-Calibration _  Linearity,  gain,  offset  and  input  bias  current 

External  Calibration  Interval _ ..... _ ..........  2  years 


Certification  and  Compliances 

CE  Mark  Compliance  c« 

for  detailed  specifications  on  environmental,  safety,  and  physical  dimensions,  please  visit  ■i.conVdigirizers. 


Figure  4.  CMRR  Versus  Frequency  for  the  PXI-5922  in  Differential  Mode 
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Holosonic  Audio  Spotlight 


audio 


^^light 


Technical 

Specifications 


Sound  Field  Distribution 


AS-16 


AS-24 


Amplifier  Specifications 


Sound  field  distribution  Is  shown 
with  equal-loudness  contours  for 
a  standard  1  kHz  tone.The  center 
area  is  loudest  at  1 00%  ampli¬ 
tude, while  the  sound  leveljust 
outside  the  Illustrated  beam  area 
Is  less  than  10%. 

Audio  Spotlight  systems  are 
much  less  sensitive  to  listener 
distance  than  traditional 
loudspeakers,  but  maximum 
performance  Is  attained  at 
roughly  1-2m  (3-6  ft)  from  the 
listener. 

Typical  levels  are  80  dB  SPL  at 
1  kHz  for  AS-16, and  85  dB  SPLfor 
AS-24  models.The  larger  AS-24 
<  1 0%  can  output  about  twice  the 

-20  dB  powerand  has  twice  low- 

frequency  range  of  the  AS-16. 

1  ft 

1  m 


Speaker  Dimensions  (thickness  ~0.5"/  1cm) 

!5  3/4'(40cm) 


Input: 

Power  draw: 

Output: 

Controls: 

Voltage: 

Dimensions: 


RCA  line-level  audio 
65W  max  (AS-24) 
50W  max  (AS-16) 

BNC  coax  cable 
(2577m  included) 
Volume,  tone,  on/off 
1 00-240V  50/60Hz 
6"wx7”dx  1.6"h 
(15cmx  1 8cm  x  4cm) 


AS-16 

New  square  models  shown. 
Legacy  round  speakers 
arealso  available. 


Holosonic  Research  Labs, Inc. 

www.holosonics.com 

6 1 7-923-4000  info@holosonics.com 


Add  sound... 

and  preserve  the  quiet. 


<•  2007  Hcfosoric  Research  Labs,  Inc.  Audio  Spodic^it  isaregstered  trademark  of  Holosorac  Research  Labs,  Inc. 


23  3/4"(«X3cm) 


Piezotronic  377B01  Condenser  Microphone 


Model  Number 
377B01 


Revision  B 
ECN  #:  25497 


PRECISION  CONDENSER  MICROPHONE 


Performance 

Nominal  Microphone  Diameter 
Frequency  Response  Characteristic  (at  0° 
incidence) 

Open  Circuit  Sensitivity  (at  250  Hz) 

Open  Circuit  Sensitivity  (+3dB)  (at  250 
Hz) 

Frequency  Range  (±2dB) 

Lower  Limiting  Frequency  (-3  dB) 
Resonant  Frequency  (90'16  hse  Shift) 
Dynamic  Range  (3%  Distortion  Limit) 
Dynamic  Range  (Cartridge  Thermal 
Noise) 

Standards  Designation  (IEC  61 094-4) 
Environmental 

Temperature  Range  (Operating) 

T emperature  Coefficient  of  Sensitivity  (14 
to  122"F,  -10  to  50’C) 

Static  Pressure  Coefficient 
Influence  of  Humidity  (0  to  95%,  non¬ 
condensing) 

I  nfluence  of  Axial  Vibration  (0. 1  g  ( 1  m/s-)) 
Electrical 

Capacitance  (Polarized) 

Polarization  Voltage 
Physical 
Housing  Material 
Venting 

Mounting  Thread  (Preamplifier) 

Mounting  Thread  (Grid) 

Size  (Diameter)  (with  grid) 

Size  (Diameter)  (without  grid) 

Size  (Height)  (with  grid) 

Size  (Height)  (without  grid) 

Weight 


4  to  80000  Hz 
0.5  to  3  Hz 
83  KHz 

>165  dB  re  20  MPa 
28  dB(A)  re  20  pPa 

WS3F 

-40  to  +248  *F 
0.004  dBTF 


60  dB  re  20  pPa 

5.1  pF 
0V 

Nickel  Alloy 
Side 

0.2244  -  60  UNS 
0.25  -  60  UNS 
0.28  in 
0.25  in 
0.41  in 
0.35  in 
0.06  oz 


4  to  80000  Hz 
0.5  to  3  Hz 
83  kHz 

>165  dB  re  20  MPa 
28  dB(A)  re  20  pPa 

WS3F 

-40  to +120  *C 
0.009  dBrc 


60  dB  re  20  pPa 

5.1  pF 
0V 

Nickel  Alloy 

5.7  mm  -  60  UNS 
6.35  mm  -  60  UNS 
7.0  mm 
6.4  mm 
10  5  mm 
9.0  mm 
1.8  gm 


All  specifications  are  at  room  temperature  unless  othemise  specified. 

In  the  interest  of  constant  product  improvement,  we  reserve  the  right  to  change  specifications  without 
notice 

ICP®  is  a  registered  trademark  of  PCB  group,  Inc. 


Optional  Versions  (Optional  versions  have  identical  specifications  and  accessories  as  listed 
for  standard  model  except  where  noted  below.  More  than  one  option  maybe  used.) 

Notes 

[1]  Typical 

[2]  Prepolarized 

Supplied  Accessories 

ACS-20  Calibration  of  Precision  Condensor  Microphones  (1) 


Entered:  LLH 

Engineer:  CCL 

Sales:  MV 

Approved:  BAM 

Spec  Number: 

Date: 

12/20/2006 

Date: 

12/21/2006 

Date: 

12/21/2006 

Date: 

12/21/2006 

31221 

® PCB  PIEZ0TR0NICS 


VIBRATION  DIVISION 


3425  Walden  Avenue 
Depew,  NY  14043 
UNITED  STATES 
Phone:  888-684-0013 
Fax:  716-685-3886 
E-mail:  vibration@pcb.com 
Web  site:  wvrw.pcb.com 


